!
! http://people.sc.fsu.edu/~jburkardt/f_src/linpack_d/linpack_d.html
!
! LINPACK
! Linear Algebra Library
! Real*8 version
!
! LINPACK is a FORTRAN90 library which solves systems of linear equations for
! a variety of matrix types and storage modes.
!
! LINPACK has officially been superseded by the LAPACK library. The LAPACK
! library uses more modern algorithms and code structure. However, the LAPACK
! library can be extraordinarily complex; what is done in a single LINPACK
! routine may correspond to 10 or 20 utility routines in LAPACK. This is fine
! if you treat LAPACK as a black box. But if you wish to learn how the
! algorithm works, or to adapt it, or to convert the code to another language,
! this is a real drawback. This is one reason I still keep a copy of LINPACK
! around.
!
! Versions of LINPACK in various arithmetic precisions are available through the
! NETLIB web site.
! Licensing:
!
! This computer code is distributed under the GNU LGPL license.
!
! List of Routines:
!
! DCHDC computes the Cholesky decomposition of a positive definite matrix.
! DCHDD downdates an augmented Cholesky decomposition.
! DCHEX updates the Cholesky factorization of a positive definite matrix.
! DCHUD updates an augmented Cholesky decomposition.
! DGBCO factors a real*8 band matrix and estimates its condition.
! DGBDI computes the determinant of a band matrix factored by DGBCO or DGBFA.
! DGBFA factors a real*8 band matrix by elimination.
! DGBSL solves a real*8 banded system factored by DGBCO or DGBFA.
! DGECO factors a real*8 matrix and estimates its condition number.
! DGEDI computes the determinant and inverse of a matrix factored by DGECO or
!       DGEFA.
! DGEFA factors a real*8 matrix.
! DGESL solves a real*8 general linear system A * X = B.
! DGTSL solves a general tridiagonal linear system.
! DPBCO factors a real*8 symmetric positive definite banded matrix.
! DPBDI computes the determinant of a matrix factored by DPBCO or DPBFA.
! DPBFA factors a real*8 symmetric positive definite matrix stored in band form.
! DPBSL solves a real*8 symmetric positive definite band factored by DPBCO or
!       DPBFA.
! DPOCO factors a real*8 symmetric positive definite matrix and estimates its
!       condition.
! DPODI computes the determinant and inverse of a certain matrix.
! DPOFA factors a real*8 symmetric positive definite matrix.
! DPOSL solves a linear system factored by DPOCO or DPOFA.
! DPPCO factors a real*8 symmetric positive definite matrix in packed form.
! DPPDI computes the determinant and inverse of a matrix factored by DPPCO or
!       DPPFA.
! DPPFA factors a real*8 symmetric positive definite matrix in packed form.
! DPPSL solves a real*8 symmetric positive definite system factored by DPPCO or
!       DPPFA.
! DPTSL solves a positive definite tridiagonal linear system.
! DQRDC computes the QR factorization of a real*8 rectangular matrix.
! DQRSL computes transformations, projections, and least squares solutions.
! DSICO factors a real*8 symmetric matrix and estimates its condition.
! DSIDI computes the determinant, inertia and inverse of a real*8 symmetric
!       matrix.
! DSIFA factors a real*8 symmetric matrix.
! DSISL solves a real*8 symmetric system factored by DSIFA.
! DSPCO factors a real*8 symmetric matrix stored in packed form.
! DSPDI computes the determinant, inertia and inverse of a real*8 symmetric
!       matrix.
! DSPFA factors a real*8 symmetric matrix stored in packed form.
! DSPSL solves the real*8 symmetric system factored by DSPFA.
! DSVDC computes the singular value decomposition of a real*8 rectangular
!       matrix.
! DTRCO estimates the condition of a real*8 triangular matrix.
! DTRDI computes the determinant and inverse of a real*8 triangular matrix.
! DTRSL solves triangular linear systems.
!
module linpack

  use mod_realkinds

  implicit none

contains

subroutine i_swap ( i, j )
!
!*******************************************************************************
!
!! I_SWAP swaps two integer values.
!
!
!  Modified:
!
!    30 November 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer I, J.  On output, the values of I and
!    J have been interchanged.
!
  implicit none
!
  integer i
  integer j
  integer k
!
  k = i
  i = j
  j = k

  return
end subroutine i_swap

function isamax ( n, x, incx )
!
!*******************************************************************************
!
!! ISAMAX finds the index of the vector element of maximum absolute value.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real(rkx) X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of SX.
!
!    Output, integer ISAMAX, the index of the element of SX of maximum
!    absolute value.
!
  implicit none
!
  integer i
  integer incx
  integer isamax
  integer ix
  integer n
  real(rkx) x(*)
  real(rkx) samax
!
  if ( n <= 0 ) then

    isamax = 0

  else if ( n == 1 ) then

    isamax = 1

  else if ( incx == 1 ) then

    isamax = 1
    samax = abs ( x(1) )

    do i = 2, n

      if ( abs ( x(i) ) > samax ) then
        isamax = i
        samax = abs ( x(i) )
      end if

    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    isamax = 1
    samax = abs ( x(ix) )

    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) > samax ) then
        isamax = i
        samax = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end function isamax

subroutine r_swap ( x, y )
!
!*******************************************************************************
!
!! R_SWAP swaps two real(rkx) values.
!
!
!  Modified:
!
!    01 May 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real(rkx) X, Y.  On output, the values of X and
!    Y have been interchanged.
!
  implicit none
!
  real(rkx) x
  real(rkx) y
  real(rkx) z
!
  z = x
  x = y
  y = z

  return
end subroutine r_swap

function samax ( n, x, incx )
!
!*******************************************************************************
!
!! SAMAX returns the maximum absolute value of the entries in a vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real(rkx) X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real(rkx) SAMAX, the maximum absolute value of an element of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer n
  real(rkx) samax
  real(rkx) x(*)
!
  if ( n <= 0 ) then

    samax = 0.0_rkx

  else if ( n == 1 ) then

    samax = abs ( x(1) )

  else if ( incx == 1 ) then

    samax = abs ( x(1) )

    do i = 2, n
      if ( abs ( x(i) ) > samax ) then
        samax = abs ( x(i) )
      end if
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    samax = abs ( x(ix) )
    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) > samax ) then
        samax = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end function samax

function sasum ( n, x, incx )
!
!*******************************************************************************
!
!! SASUM sums the absolute values of the entries of a vector.
!
!
!  Modified:
!
!    15 February 2001
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real(rkx) X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!    INCX must not be negative.
!
!    Output, real(rkx) SASUM, the sum of the absolute values of X.
!
  implicit none
!
  integer incx
  integer n
  real(rkx) sasum
  real(rkx) x(*)
!
  sasum = sum ( abs ( x(1:1+(n-1)*incx:incx) ) )

  return
end function sasum

subroutine saxpy ( n, sa, x, incx, y, incy )
!
!*******************************************************************************
!
!! SAXPY adds a constant times one vector to another.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real(rkx) SA, the multiplier.
!
!    Input, real(rkx) X(*), the vector to be scaled and added to Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real(rkx) Y(*), the vector to which a multiple of X is to
!    be added.
!
!    Input, integer INCY, the increment between successive entries of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real(rkx) sa
  real(rkx) x(*)
  real(rkx) y(*)
!
  if ( n <= 0 ) then

  else if ( sa == 0.0_rkx ) then

  else if ( incx == 1 .and. incy == 1 ) then

    y(1:n) = y(1:n) + sa * x(1:n)

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      y(iy) = y(iy) + sa * x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end subroutine saxpy

subroutine schdc ( a, lda, p, work, ipvt, job, info )
!
!*******************************************************************************
!
!! SCHDC computes the Cholesky decomposition of a positive definite matrix.
!
!
!  Discussion:
!
!    A pivoting option allows the user to estimate the condition of a
!    positive definite matrix or determine the rank of a positive
!    semidefinite matrix.
!
!    For positive definite matrices, INFO = P is the normal return.
!
!    For pivoting with positive semidefinite matrices, INFO will
!    in general be less than P.  However, INFO may be greater than
!    the rank of A, since rounding error can cause an otherwise zero
!    element to be positive.  Indefinite systems will always cause
!    INFO to be less than P.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Author:
!
!    J Dongarra and G Stewart,
!    Argonne National Laboratory and University of Maryland.
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,P).
!    On input, A contains the matrix whose decomposition is to
!    be computed.  Only the upper half of A need be stored.
!    The lower part of the array a is not referenced.
!    On output, A contains in its upper half the Cholesky factor
!    of the input matrix, as it has been permuted by pivoting.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer P, the order of the matrix.
!
!    Input, real(rkx) WORK(P) is a work array.
!
!    Input/output, integer IPVT(P).
!    On input, IPVT contains integers that control the selection
!    of the pivot elements, if pivoting has been requested.
!    Each diagonal element A(K,K) is placed in one of three classes
!    according to the value of IPVT(K).
!
!      > 0, then X(K) is an initial element.
!      = 0, then X(K) is a free element.
!      < 0, then X(K) is a final element.
!
!    Before the decomposition is computed, initial elements are moved by
!    symmetric row and column interchanges to the beginning of the array A
!    and final elements to the end.  Both initial and final elements are
!    frozen in place during the computation and only free elements are moved.
!    At the K-th stage of the reduction, if A(K,K) is occupied by a free
!    element, it is interchanged with the largest free element A(L,L) with
!    L >= K.  IPVT is not referenced if JOB is 0.
!
!    On output, IPVT(J) contains the index of the diagonal element
!    of A that was moved into the J-th position, if pivoting was requested.
!
!    Input, integer JOB, initiates column pivoting.
!    0, no pivoting is done.
!    nonzero, pivoting is done.
!
!    Output, integer INFO, contains the index of the last positive diagonal
!    element of the Cholesky factor.
!
  implicit none
!
  integer lda
  integer p
!
  real(rkx) a(lda,p)
  integer info
  integer j
  integer job
  integer ipvt(p)
  integer k
  integer l
  real(rkx) maxdia
  integer maxl
  logical negk
  integer pl
  integer pu
  logical swapk
  real(rkx) temp
  real(rkx) work(*)
!
  pl = 1
  pu = 0
  info = p

  if ( job /= 0 ) then
!
!  Pivoting has been requested.
!  Rearrange the the elements according to IPVT.
!
    do k = 1, p

      swapk = ipvt(k) > 0

      negk = ipvt(k) < 0

      if ( negk ) then
        ipvt(k) = -k
      else
        ipvt(k) = k
      end if

      if ( swapk ) then

        if ( k /= pl ) then

          call sswap ( pl-1, a(1,k), 1, a(1,pl), 1 )

          call r_swap ( a(k,k), a(pl,pl) )

          do j = pl+1, p

            if ( j < k ) then
              call r_swap ( a(pl,j), a(j,k) )
            else if ( j > k ) then
              call r_swap ( a(pl,j), a(k,j) )
            end if

          end do

          ipvt(k) = ipvt(pl)
          ipvt(pl) = k

        end if

        pl = pl + 1

      end if

    end do

    pu = p

    do k = p, pl, -1

      if ( ipvt(k) < 0 ) then

        ipvt(k) = -ipvt(k)

        if ( pu /= k ) then

          call sswap ( k-1, a(1,k), 1, a(1,pu), 1 )
          call r_swap ( a(k,k), a(pu,pu) )

          do j = k+1, p

            if ( j < pu ) then
              call r_swap ( a(k,j), a(j,pu) )
            else if ( j > pu ) then
              call r_swap ( a(k,j), a(pu,j) )
            end if

          end do

          call i_swap ( ipvt(k), ipvt(pu) )

        end if

        pu = pu - 1

      end if

    end do

  end if

  do k = 1, p
!
!  Reduction loop.
!
    maxdia = a(k,k)
    maxl = k
!
!  Determine the pivot element.
!
    if ( k >= pl .and. k < pu ) then

      do l = k+1, pu
        if ( a(l,l) > maxdia ) then
          maxdia = a(l,l)
          maxl = l
        end if
      end do

    end if
!
!  Quit if the pivot element is not positive.
!
    if ( maxdia <= 0.0_rkx ) then
      info = k - 1
      return
    end if
!
!  Start the pivoting and update IPVT.
!
    if ( k /= maxl ) then

      call sswap ( k-1, a(1,k), 1, a(1,maxl), 1 )

      a(maxl,maxl) = a(k,k)
      a(k,k) = maxdia

      call i_swap ( ipvt(maxl), ipvt(k) )

    end if
!
!  Reduction step.
!  Pivoting is contained across the rows.
!
    work(k) = sqrt ( a(k,k) )
    a(k,k) = work(k)

    do j = k+1, p

      if ( k /= maxl ) then

        if ( j < maxl ) then
          call r_swap ( a(k,j), a(j,maxl) )
        else if ( j > maxl ) then
          call r_swap ( a(k,j), a(maxl,j) )
        end if

      end if

      a(k,j) = a(k,j) / work(k)
      work(j) = a(k,j)
      temp = -a(k,j)
      call saxpy ( j-k, temp, work(k+1), 1, a(k+1,j), 1 )

    end do

  end do

  return
end subroutine schdc

subroutine schdd ( r, ldr, p, x, z, ldz, nz, y, rho, c, s, info )
!
!*******************************************************************************
!
!! SCHDD downdates an augmented Cholesky decomposition.
!
!
!  Discussion:
!
!    SCHDD can also downdate the triangular factor of an augmented QR
!    decomposition.
!
!    Specifically, given an upper triangular matrix R of order P, a
!    row vector X, a column vector Z, and a scalar Y, SCHDD
!    determines an orthogonal matrix U and a scalar ZETA such that
!
!          ( R   Z  )     ( RR  ZZ )
!      U * (        )  =  (        ),
!          ( 0 ZETA )     (  X   Y )
!
!    where RR is upper triangular.
!
!    If R and Z have been obtained from the factorization of a least squares
!    problem, then RR and ZZ are the factors corresponding to the problem
!    with the observation (X,Y) removed.  In this case, if RHO
!    is the norm of the residual vector, then the norm of
!    the residual vector of the downdated problem is
!    sqrt ( RHO**2 - ZETA**2 ). SCHDD will simultaneously downdate
!    several triplets (Z, Y, RHO) along with R.
!
!    For a less terse description of what SCHDD does and how
!    it may be applied, see the LINPACK guide.
!
!    The matrix U is determined as the product U(1)*...*U(P)
!    where U(I) is a rotation in the (P+1,I)-plane of the form
!
!      ( C(I)      -S(I)    )
!      (                    ).
!      ( S(I)       C(I)    )
!
!    The rotations are chosen so that C(I) is real.
!
!    The user is warned that a given downdating problem may be impossible
!    to accomplish or may produce inaccurate results.  For example, this
!    can happen if X is near a vector whose removal will reduce the
!    rank of R.  Beware.
!
!  Modified:
!
!    17 April 2002
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) R(LDR,P), the upper triangular matrix that is to be
!    downdated.  The part of R below the diagonal is not referenced.
!
!    Input, integer LDR, the leading dimension of the array R.
!    LDR must be at least P.
!
!    Input, integer P, the order of the matrix R.
!
!    Input, real(rkx) X(P), the row vector that is to be removed from R.
!
!    Input/output, real(rkx) Z(LDZ,NZ), an array of NZ P-vectors which are to
!    be downdated along with R.
!
!    Input, integer LDZ, the leading dimension of the array Z.
!    LDZ must be at least P.
!
!    Input, integer NZ, the number of vectors to be downdated.
!    NZ may be zero, in which case Z, Y, and RHO are not referenced.
!
!    Input, real(rkx) Y(NZ), the scalars for the downdating of the vectors Z.
!
!    Input/output, real(rkx) RHO(NZ), the norms of the residual vectors.  On
!    output these have been changed along with R and Z.
!
!    Output, real(rkx) C(P), S(P), the cosines and sines of the transforming
!    rotations.
!
!    Output, integer INFO, return flag.
!     0, the entire downdating was successful.
!    -1, if R could not be downdated.  In this case, all quantities
!        are left unaltered.
!     1, if some RHO could not be downdated.  The offending RHO's are
!        set to -1.
!
  implicit none
!
  integer ldr
  integer ldz
  integer nz
  integer p
!
  real(rkx) a
  real(rkx) alpha
  real(rkx) azeta
  real(rkx) b
  real(rkx) c(p)
  integer i
  integer ii
  integer info
  integer j
  real(rkx) norm
  real(rkx) r(ldr,p)
  real(rkx) rho(nz)
  real(rkx) s(p)
  real(rkx) scale
  real(rkx) t
  real(rkx) x(p)
  real(rkx) xx
  real(rkx) y(nz)
  real(rkx) z(ldz,nz)
  real(rkx) zeta
!
!  Solve the system R'*A = X, placing the result in the array S.
!
  info = 0
  s(1) = x(1) / r(1,1)

  do j = 2, p
    s(j) = x(j) - dot_product ( r(1:j-1,j), s(1:j-1) )
    s(j) = s(j) / r(j,j)
  end do

  norm = snrm2 ( p, s, 1 )

  if ( norm >= 1.0_rkx ) then
    info = -1
    return
  end if

  alpha = sqrt ( 1.0_rkx - norm**2 )
!
!  Determine the transformations.
!
  do ii = 1, p
    i = p - ii + 1
    scale = alpha + abs ( s(i) )
    a = alpha / scale
    b = s(i) / scale
    norm = sqrt ( a**2 + b**2 )
    c(i) = a / norm
    s(i) = b / norm
    alpha = scale * norm
  end do
!
!  Apply the transformations to R.
!
  do j = 1, p
    xx = 0.0_rkx
    do ii = 1, j
      i = j - ii + 1
      t = c(i) * xx + s(i) * r(i,j)
      r(i,j) = c(i) * r(i,j) - s(i) * xx
      xx = t
    end do
  end do
!
!  If required, downdate Z and RHO.
!
  do j = 1, nz

    zeta = y(j)
    do i = 1, p
      z(i,j) = ( z(i,j) - s(i) * zeta ) / c(i)
      zeta = c(i) * zeta - s(i) * z(i,j)
    end do

    azeta = abs ( zeta )

    if ( azeta > rho(j) ) then
      info = 1
      rho(j) = -1.0_rkx
    else
      rho(j) = rho(j) * sqrt ( 1.0_rkx - ( azeta / rho(j) )**2 )
    end if

  end do

  return
end subroutine schdd

subroutine schex ( r, ldr, p, k, l, z, ldz, nz, c, s, job )
!
!*******************************************************************************
!
!! SCHEX updates the Cholesky factorization of a positive definite matrix.
!
!
!  Discussion:
!
!    The factorization has the form
!
!      A = R' * R
!
!    where A is a positive definite matrix of order P.
!
!    The updating involves diagonal permutations of the form
!
!      E' * A * E
!
!    where E is a permutation matrix.  Specifically, given
!    an upper triangular matrix R and a permutation matrix
!    E (which is specified by K, L, and JOB), SCHEX determines
!    a orthogonal matrix U such that
!
!      U * R * E = RR,
!
!    where RR is upper triangular.  At the user's option, the
!    transformation U will be multiplied into the array Z.
!    If A = X'*X, so that R is the triangular part of the
!    QR factorization of X, then RR is the triangular part of the
!    QR factorization of X*E, i.e. X with its columns permuted.
!    For a less terse description of what SCHEX does and how
!    it may be applied, see the LINPACK guide.
!
!    The matrix Q is determined as the product U(L-K)*...*U(1)
!    of plane rotations of the form
!
!      (    C(I)       S(I) )
!      (                    ),
!      (   -S(I)       C(I) )
!
!    where C(I) is real, the rows these rotations operate on
!    are described below.
!
!    There are two types of permutations, which are determined
!    by the value of JOB.
!
!    1, right circular shift.  The columns are rearranged in the order:
!
!         1,...,K-1,L,K,K+1,...,L-1,L+1,...,P.
!
!       U is the product of L-K rotations U(I), where U(I)
!       acts in the (L-I,L-I+1)-plane.
!
!    2, left circular shift: the columns are rearranged in the order
!
!         1,...,K-1,K+1,K+2,...,L,K,L+1,...,P.
!
!       U is the product of L-K rotations U(I), where U(I)
!       acts in the (K+I-1,K+I)-plane.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) R(LDR,P).  On input, the upper triangular factor
!    that is to be updated.  Elements of R below the diagonal are not
!    referenced.  On output, R has been updated.
!
!    Input, integer LDR, the leading dimension of the array R.
!    LDR must be at least P.
!
!    Input, integer P, the order of the matrix R.
!
!    Input, integer K, the first column to be permuted.
!
!    Input, integer L, the last column to be permuted.
!    L must be strictly greater than K.
!
!    Input/output real(rkx) Z(LDZ,NZ), an array of NZ P-vectors into which the
!    transformation U is multiplied.  Z is not referenced if NZ = 0.
!    On output, Z has been updated.
!
!    Input, integer LDZ, the leading dimension of the array Z.
!    LDZ must be at least P.
!
!    Input, integer NZ, the number of columns of the matrix Z.
!
!    Input, integer JOB, determines the type of permutation.
!    1, right circular shift.
!    2, left circular shift.
!
!    Output, real(rkx) C(P), S(P), the cosines and sines of the transforming
!    rotations.
!
  implicit none
!
  integer ldr
  integer ldz
  integer p
  integer nz
!
  real(rkx) c(p)
  integer job
  integer i
  integer ii
  integer il
  integer iu
  integer j
  integer jj
  integer k
  integer l
  integer lm1
  integer lmk
  real(rkx) r(ldr,p)
  real(rkx) s(p)
  real(rkx) t
  real(rkx) z(ldz,nz)
!
!  Initialize
!
  lmk = l - k
  lm1 = l - 1
!
!  Right circular shift.
!
  if ( job == 1 ) then
!
!  Reorder the columns.
!
     do i = 1, l
       ii = l - i + 1
       s(i) = r(ii,l)
     end do

     do jj = k, lm1
       j = lm1 - jj + k
       do i = 1, j
         r(i,j+1) = r(i,j)
       end do
       r(j+1,j+1) = 0.0_rkx
     end do

     do i = 1, k-1
       ii = l - i + 1
       r(i,k) = s(ii)
     end do
!
!  Calculate the rotations.
!
     t = s(1)
     do i = 1, lmk
       call srotg ( s(i+1), t, c(i), s(i) )
       t = s(i+1)
     end do

     r(k,k) = t

     do j = k+1, p
       il = max (1,l-j+1)
       do ii = il, lmk
         i = l - ii
         t = c(ii) * r(i,j) + s(ii) * r(i+1,j)
         r(i+1,j) = c(ii) * r(i+1,j) - s(ii) * r(i,j)
         r(i,j) = t
       end do
     end do
!
!  If required, apply the transformations to Z.
!
     do j = 1, nz
       do ii = 1, lmk
         i = l - ii
         t = c(ii) * z(i,j) + s(ii) * z(i+1,j)
         z(i+1,j) = c(ii) * z(i+1,j) - s(ii) * z(i,j)
         z(i,j) = t
       end do
     end do
!
!  Left circular shift.
!
  else
!
!  Reorder the columns.
!
     do i = 1, k
       ii = lmk + i
       s(ii) = r(i,k)
     end do

     do j = k, lm1
       do i = 1, j
         r(i,j) = r(i,j+1)
       end do
       jj = j - k + 1
       s(jj) = r(j+1,j+1)
     end do

     do i = 1, k
       ii = lmk + i
       r(i,l) = s(ii)
     end do

     r(k+1:l,l) = 0.0_rkx
!
!  Reduction loop.
!
     do j = k, p
!
!  Apply the rotations.
!
        if ( j /= k ) then

           iu = min ( j-1, l-1 )

           do i = k, iu
              ii = i - k + 1
              t = c(ii) * r(i,j) + s(ii) * r(i+1,j)
              r(i+1,j) = c(ii) * r(i+1,j) - s(ii) * r(i,j)
              r(i,j) = t
           end do

        end if

        if ( j < l ) then
           jj = j - k + 1
           t = s(jj)
           call srotg ( r(j,j), t, c(jj), s(jj) )
        end if

     end do
!
!  Apply the rotations to Z.
!
     do j = 1, nz
        do i = k, lm1
           ii = i - k + 1
           t = c(ii) * z(i,j) + s(ii) * z(i+1,j)
           z(i+1,j) = c(ii) * z(i+1,j) - s(ii) * z(i,j)
           z(i,j) = t
        end do
     end do

  end if

  return
end subroutine schex

subroutine schud ( r, ldr, p, x, z, ldz, nz, y, rho, c, s )
!
!*******************************************************************************
!
!! SCHUD updates an augmented Cholesky decomposition.
!
!
!  Discussion:
!
!    SCHUD can also update the triangular part of an augmented QR
!    decomposition.
!
!    Specifically, given an upper triangular matrix R of order P, a row vector
!    X, a column vector Z, and a scalar Y, SCHUD determines a unitary matrix
!    U and a scalar ZETA such that
!
!           (R  Z)     (RR   ZZ )
!      U  * (    )  =  (        ),
!           (X  Y)     ( 0  ZETA)
!
!    where RR is upper triangular.
!
!    If R and Z have been obtained from the factorization of a least squares
!    problem, then RR and ZZ are the factors corresponding to the problem
!    with the observation (X,Y) appended.  In this case, if RHO is the
!    norm of the residual vector, then the norm of the residual vector of
!    the updated problem is sqrt ( RHO**2 + ZETA**2 ).  SCHUD will
!    simultaneously update several triplets (Z, Y, RHO).
!
!    For a less terse description of what SCHUD does and how
!    it may be applied, see the LINPACK guide.
!
!    The matrix U is determined as the product U(P)*...*U(1),
!    where U(I) is a rotation in the (I,P+1) plane of the form
!
!      (     C(I)      S(I) )
!      (                    ).
!      (    -S(I)      C(I) )
!
!    The rotations are chosen so that C(I) is real.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) R(LDR,P), the upper triangular matrix that is to be updated.
!    the part of R below the diagonal is not referenced.
!
!    Input, integer LDR, the leading dimension of the array R.
!    LDR must be at least equal to P.
!
!    Input, integer P, the order of the matrix R.
!
!    Input, real(rkx) X(P), the row to be added to R.
!
!    Input/output, real(rkx) Z(LDZ,NZ), contains NZ P-vectors to be updated with R.
!
!    Input, integer LDZ, the leading dimension of the array Z.
!    LDZ must be at least P.
!
!    Input, integer NZ, the number of vectors to be updated.  NZ may be
!    zero, in which case Z, Y, and RHO are not referenced.
!
!    Input, real(rkx) Y(NZ), the scalars for updating the vectors Z.
!
!    Input/output, real(rkx) RHO(NZ).  On input, the norms of the residual
!    vectors that are to be updated.  If RHO(J) is negative, it is left
!    unaltered.
!
!    Output, real(rkx) C(P), S(P), the cosines and sines of the transforming
!    rotations.
!
  implicit none
!
  integer ldr
  integer ldz
  integer nz
  integer p
!
  real(rkx) azeta
  real(rkx) c(p)
  integer i
  integer j
  real(rkx) r(ldr,p)
  real(rkx) rho(nz)
  real(rkx) s(p)
  real(rkx) scale
  real(rkx) t
  real(rkx) x(p)
  real(rkx) xj
  real(rkx) y(nz)
  real(rkx) z(ldz,nz)
  real(rkx) zeta
!
!  Update R.
!
  do j = 1, p

     xj = x(j)
!
!  Apply the previous rotations.
!
     do i = 1, j-1
        t = c(i) * r(i,j) + s(i) * xj
        xj = c(i) * xj - s(i) * r(i,j)
        r(i,j) = t
     end do
!
!  Compute the next rotation.
!
     call srotg ( r(j,j), xj, c(j), s(j) )

  end do
!
!  If required, update Z and RHO.
!
  do j = 1, nz

     zeta = y(j)

     do i = 1, p
        t = c(i) * z(i,j) + s(i) * zeta
        zeta = c(i) * zeta - s(i) * z(i,j)
        z(i,j) = t
     end do

     azeta = abs ( zeta )

     if ( azeta /= 0.0_rkx .and. rho(j) >= 0.0 ) then
        scale = azeta + rho(j)
        rho(j) = scale * sqrt ( ( azeta / scale )**2 + ( rho(j) / scale )**2 )
     end if

  end do

  return
end subroutine schud

subroutine scopy ( n, x, incx, y, incy )
!
!*******************************************************************************
!
!! SCOPY copies one real(rkx) vector into another.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real(rkx) X(*), the vector to be copied into Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real(rkx) Y(*), the copy of X.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real(rkx) x(*)
  real(rkx) y(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    y(1:n) = x(1:n)

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      y(iy) = x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end subroutine scopy

function sdot ( n, x, incx, y, incy )
!
!*******************************************************************************
!
!! SDOT forms the dot product of two vectors.
!
!
!  Discussion:
!
!    For many purposes, SDOT can be replaced by the FORTRAN 90
!    intrinsic DOT_PRODUCT.
!
!  Modified:
!
!    02 June 2000
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input, real(rkx) X(*), one of the vectors to be multiplied.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input, real(rkx) Y(*), one of the vectors to be multiplied.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
!    Output, real(rkx) SDOT, the dot product of X and Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real(rkx) sdot
  real(rkx) stemp
  real(rkx) x(*)
  real(rkx) y(*)
!
  if ( n <= 0 ) then

    sdot = 0.0_rkx

  else if ( incx == 1 .and. incy == 1 ) then

    sdot = dot_product ( x(1:n), y(1:n) )

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    stemp = 0.0_rkx
    do i = 1, n
      stemp = stemp + x(ix) * y(iy)
      ix = ix + incx
      iy = iy + incy
    end do

    sdot = stemp

  end if

  return
end function sdot

subroutine sgbco ( abd, lda, n, ml, mu, ipvt, rcond, z )
!
!*******************************************************************************
!
!! SGBCO factors a real(rkx) band matrix and estimates its condition.
!
!
!  Discussion:
!
!    If RCOND is not needed, SGBFA is slightly faster.
!
!    To solve A*X = B, follow SGBCO by SGBSL.
!
!    To compute inverse(A)*C, follow SGBCO by SGBSL.
!
!    To compute determinant(A), follow SGBCO by SGBDI.
!
!  Example:
!
!    If the original matrix is
!
!      11 12 13  0  0  0
!      21 22 23 24  0  0
!       0 32 33 34 35  0
!       0  0 43 44 45 46
!       0  0  0 54 55 56
!       0  0  0  0 65 66
!
!    then for proper band storage,
!
!      N = 6, ML = 1, MU = 2, LDA >= 5  and ABD should contain
!
!       *  *  *  +  +  +      * = not used
!       *  * 13 24 35 46      + = used for pivoting
!       * 12 23 34 45 56
!      11 22 33 44 55 66
!      21 32 43 54 65  *
!
!  Band storage:
!
!    If A is a band matrix, the following program segment
!    will set up the input.
!
!      ml = (band width below the diagonal)
!      mu = (band width above the diagonal)
!      m = ml + mu + 1
!
!      do j = 1, n
!        i1 = max ( 1, j-mu )
!        i2 = min ( n, j+ml )
!        do i = i1, i2
!          k = i - j + m
!          abd(k,j) = a(i,j)
!        end do
!      end do
!
!    This uses rows ML+1 through 2*ML+MU+1 of ABD.  In addition, the first
!    ML rows in ABD are used for elements generated during the
!    triangularization.
!
!    The total number of rows needed in ABD is  2*ML+MU+1.  The ML+MU by
!    ML+MU upper left triangle and the ML by ML lower right triangle are
!    not referenced.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) ABD(LDA,N).  On input, the matrix in band storage.
!    The columns of the matrix are stored in the columns of ABD and
!    the diagonals of the matrix are stored in rows ML+1 through 2*ML+MU+1
!    of ABD.  On output, an upper triangular matrix in band storage and
!    the multipliers which were used to obtain it.  The factorization can
!    be written A = L*U where L is a product of permutation and unit lower
!    triangular matrices and U is upper triangular.
!
!    Input, integer LDA, the leading dimension of the array ABD.
!    LDA must be >= 2*ML + MU + 1.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer ML, MU, the number of diagonals below and above the
!    main diagonal.  0 <= ML < N, 0 <= MU < N.
!
!    Output, integer IPVT(N), the pivot indices.
!
!    Output, real(rkx) RCOND, an estimate of the reciprocal condition of A.
!    For the system A*X = B, relative perturbations in A and B of size
!    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
!    If RCOND is so small that the logical expression
!      1.0_rkx + RCOND == 1.0_rkx
!    is true, then A may be singular to working precision.  In particular,
!    RCOND is zero if exact singularity is detected or the estimate underflows.
!
!    Workspace, real(rkx) Z(N), a work vector whose contents are usually unimportant.
!    If A is close to a singular matrix, then Z is an approximate null vector
!    in the sense that
!      norm(A*Z) = RCOND * norm(A) * norm(Z).
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) abd(lda,n)
  real(rkx) anorm
  real(rkx) ek
  integer info
  integer ipvt(n)
  integer is
  integer j
  integer ju
  integer k
  integer l
  integer la
  integer lm
  integer lz
  integer m
  integer ml
  integer mm
  integer mu
  real(rkx) rcond
  real(rkx) s
  real(rkx) sm
  real(rkx) t
  real(rkx) wk
  real(rkx) wkm
  real(rkx) ynorm
  real(rkx) z(n)
!
!  Compute the 1-norm of A.
!
  anorm = 0.0_rkx
  l = ml + 1
  is = l + mu
  do j = 1, n
    anorm = max ( anorm, sum ( abs ( abd(is:is+l-1,j) ) ) )
    if ( is > ml + 1 ) is = is - 1
    if ( j <= mu ) l = l + 1
    if ( j >= n - ml ) l = l - 1
  end do
!
!  Factor.
!
  call sgbfa ( abd, lda, n, ml, mu, ipvt, info )
!
!  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
!
!  Estimate = norm(Z)/norm(Y) where  a*z = y  and A'*Y = E.
!
!  A' is the transpose of A.  The components of E are
!  chosen to cause maximum local growth in the elements of W where
!  U'*W = E.  The vectors are frequently rescaled to avoid
!  overflow.
!
!  Solve U'*W = E.
!
  ek = 1.0_rkx
  z(1:n) = 0.0_rkx
  m = ml + mu + 1
  ju = 0

  do k = 1, n

     if ( z(k) /= 0.0_rkx ) then
       ek = sign ( ek, -z(k) )
     end if

     if ( abs ( ek - z(k) ) > abs ( abd(m,k) ) ) then
       s = abs ( abd(m,k) ) / abs ( ek - z(k) )
       z(1:n) = s * z(1:n)
       ek = s * ek
     end if

     wk = ek - z(k)
     wkm = -ek - z(k)
     s = abs ( wk )
     sm = abs ( wkm )

     if ( abd(m,k) /= 0.0_rkx ) then
       wk = wk / abd(m,k)
       wkm = wkm / abd(m,k)
     else
       wk = 1.0_rkx
       wkm = 1.0_rkx
     end if

     ju = min ( max ( ju, mu+ipvt(k) ), n )
     mm = m

     if ( k+1 <= ju ) then

        do j = k+1, ju
          mm = mm - 1
          sm = sm + abs ( z(j ) + wkm * abd(mm,j) )
          z(j) = z(j) + wk * abd(mm,j)
          s = s + abs ( z(j) )
        end do

        if ( s < sm ) then
           t = wkm - wk
           wk = wkm
           mm = m
           do j = k+1, ju
             mm = mm - 1
             z(j) = z(j) + t * abd(mm,j)
           end do
        end if

     end if

     z(k) = wk

  end do

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )
!
!  Solve L'*Y = W.
!
  do k = n, 1, -1

     lm = min ( ml, n-k )

     if ( k < n ) then
       z(k) = z(k) + dot_product ( abd(m+1:m+lm,k), z(k+1:k+lm) )
     end if

     if ( abs ( z(k) ) > 1.0_rkx ) then
       s = 1.0_rkx / abs ( z(k) )
       z(1:n) = s * z(1:n)
     end if

     l = ipvt(k)
     t = z(l)
     z(l) = z(k)
     z(k) = t
  end do

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )
  ynorm = 1.0_rkx
!
!  Solve L*V = Y.
!
  do k = 1, n

     l = ipvt(k)
     t = z(l)
     z(l) = z(k)
     z(k) = t
     lm = min ( ml, n-k )

     if ( k < n ) then
       call saxpy ( lm, t, abd(m+1,k), 1, z(k+1), 1 )
     end if

     if ( abs ( z(k) ) > 1.0_rkx ) then
       s = 1.0_rkx / abs ( z(k) )
       z(1:n) = s * z(1:n)
       ynorm = s * ynorm
     end if

  end do

  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = s * ynorm
!
!  Solve U*Z = W.
!
  do k = n, 1, -1

     if ( abs ( z(k) ) > abs ( abd(m,k) ) ) then
       s = abs ( abd(m,k) ) / abs ( z(k) )
       z(1:n) = s * z(1:n)
       ynorm = s * ynorm
     end if

     if ( abd(m,k) /= 0.0_rkx ) then
       z(k) = z(k) / abd(m,k)
     else
       z(k) = 1.0_rkx
     end if

     lm = min ( k, m ) - 1
     la = m - lm
     lz = k - lm
     t = -z(k)
     call saxpy ( lm, t, abd(la,k), 1, z(lz), 1 )

  end do
!
!  Make ZNORM = 1.0.
!
  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = s * ynorm

  if ( anorm /= 0.0_rkx ) then
    rcond = ynorm / anorm
  else
    rcond = 0.0_rkx
  end if

  return
end subroutine sgbco

subroutine sgbdi ( abd, lda, n, ml, mu, ipvt, det )
!
!*******************************************************************************
!
!! SGBDI computes the determinant of a band matrix factored by SGBCO or SGBFA.
!
!
!  Discussion:
!
!    If the inverse is needed, use SGBSL N times.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) ABD(LDA,N), the output from SGBCO or SGBFA.
!
!    Input, integer LDA, the leading dimension of the array ABD.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer ML, MU, the number of diagonals below and above the
!    main diagonal.  0 <= ML < N, 0 <= MU < N.
!
!    Input, integer IPVT(N), the pivot vector from SGBCO or SGBFA.
!
!    Output, real(rkx) DET(2), the determinant of the original matrix.
!      determinant = DET(1) * 10.0**DET(2)
!    with  1.0_rkx <= abs ( DET(1) ) < 10.0_rkx or DET(1) = 0.0_rkx.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) abd(lda,n)
  real(rkx) det(2)
  integer i
  integer ipvt(n)
  integer m
  integer ml
  integer mu
  real(rkx), parameter :: ten = 10.0_rkx
!
  m = ml + mu + 1
  det(1) = 1.0_rkx
  det(2) = 0.0_rkx

  do i = 1, n

     if ( ipvt(i) /= i ) then
       det(1) = -det(1)
     end if

     det(1) = abd(m,i) * det(1)

     if ( det(1) == 0.0_rkx ) then
       return
     end if

     do while ( abs ( det(1) ) < 1.0_rkx )
       det(1) = ten * det(1)
       det(2) = det(2) - 1.0_rkx
     end do

     do while ( abs ( det(1) ) >= ten )
       det(1) = det(1) / ten
       det(2) = det(2) + 1.0_rkx
     end do

  end do

  return
end subroutine sgbdi

subroutine sgbfa ( abd, lda, n, ml, mu, ipvt, info )
!
!*******************************************************************************
!
!! SGBFA factors a real(rkx) band matrix by elimination.
!
!
!  Discussion:
!
!    SGBFA is usually called by SGBCO, but it can be called
!    directly with a saving in time if RCOND is not needed.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) ABD(LDA,N).  On input, contains the matrix in band
!    storage.  The columns of the matrix are stored in the columns of ABD
!    and the diagonals of the matrix are stored in rows ML+1 through
!    2*ML+MU+1 of ABD.  On output, an upper triangular matrix in band storage
!    and the multipliers which were used to obtain it.  The factorization
!    can be written A = L*U where L is a product of permutation and unit lower
!    triangular matrices and U is upper triangular.
!
!    Input, integer LDA, the leading dimension of the array ABD.
!    LDA must be >= 2*ML + MU + 1.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer ML, MU, the number of diagonals below and above the
!    main diagonal.  0 <= ML < N, 0 <= MU < N.
!
!    Output, integer IPVT(N), the pivot indices.
!
!    Output, integer INFO, error flag.
!    0, normal value.
!    K, if U(K,K) == 0.0_rkx.  This is not an error condition for this
!      subroutine, but it does indicate that SGBSL will divide by zero if
!      called.  Use RCOND in SGBCO for a reliable indication of singularity.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) abd(lda,n)
  integer i0
  integer info
  integer ipvt(n)
  integer j
  integer j0
  integer j1
  integer ju
  integer jz
  integer k
  integer l
  integer lm
  integer m
  integer ml
  integer mm
  integer mu
  real(rkx) t
!
  m = ml + mu + 1
  info = 0
!
!  Zero initial fill-in columns.
!
  j0 = mu + 2
  j1 = min ( n, m ) - 1

  do jz = j0, j1
    i0 = m + 1 - jz
    abd(i0:ml,jz) = 0.0_rkx
  end do

  jz = j1
  ju = 0
!
!  Gaussian elimination with partial pivoting.
!
  do k = 1, n-1
!
!  Zero next fill-in column.
!
     jz = jz + 1
     if ( jz <= n ) then
       abd(1:ml,jz) = 0.0_rkx
     end if
!
!  Find L = pivot index.
!
     lm = min ( ml, n-k )
     l = isamax ( lm+1, abd(m,k), 1 ) + m - 1
     ipvt(k) = l + k - m
!
!  Zero pivot implies this column already triangularized.
!
     if ( abd(l,k) == 0.0_rkx ) then

       info = k
!
!  Interchange if necessary.
!
     else

        if ( l /= m ) then
          call r_swap ( abd(l,k), abd(m,k) )
        end if
!
!  Compute multipliers.
!
        t = -1.0_rkx / abd(m,k)
        call sscal ( lm, t, abd(m+1,k), 1 )
!
!  Row elimination with column indexing.
!
        ju = min ( max ( ju, mu+ipvt(k) ), n )
        mm = m

        do j = k+1, ju
           l = l - 1
           mm = mm - 1
           t = abd(l,j)
           if ( l /= mm ) then
              abd(l,j) = abd(mm,j)
              abd(mm,j) = t
           end if
           call saxpy ( lm, t, abd(m+1,k), 1, abd(mm+1,j), 1 )
        end do

     end if

  end do

  ipvt(n) = n

  if ( abd(m,n) == 0.0_rkx ) then
    info = n
  end if

  return
end subroutine sgbfa

subroutine sgbsl ( abd, lda, n, ml, mu, ipvt, b, job )
!
!*******************************************************************************
!
!! SGBSL solves a real(rkx) banded system factored by SGBCO or SGBFA.
!
!
!  Discussion:
!
!    SGBSL can solve either A * X = B  or  A' * X = B.
!
!    A division by zero will occur if the input factor contains a
!    zero on the diagonal.  Technically this indicates singularity
!    but it is often caused by improper arguments or improper
!    setting of LDA.  It will not occur if the subroutines are
!    called correctly and if SGBCO has set RCOND > 0.0_rkx
!    or SGBFA has set INFO == 0.
!
!    To compute inverse(A) * C  where C is a matrix with P columns:
!
!      call sgbco ( abd, lda, n, ml, mu, ipvt, rcond, z )
!
!      if ( rcond is too small ) then
!        exit
!      end if
!
!      do j = 1, p
!        call sgbsl ( abd, lda, n, ml, mu, ipvt, c(1,j), 0 )
!      end do
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) ABD(LDA,N), the output from SGBCO or SGBFA.
!
!    Input, integer LDA, the leading dimension of the array ABD.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer ML, MU, the number of diagonals below and above the
!    main diagonal.  0 <= ML < N, 0 <= MU < N.
!
!    Input, integer IPVT(N), the pivot vector from SGBCO or SGBFA.
!
!    Input/output, real(rkx) B(N).  On input, the right hand side.
!    On output, the solution.
!
!    Input, integer JOB, job choice.
!    0, solve A*X=B.
!    nonzero, solve A'*X=B.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) abd(lda,n)
  real(rkx) b(n)
  integer ipvt(n)
  integer job
  integer k
  integer l
  integer la
  integer lb
  integer lm
  integer m
  integer ml
  integer mu
  real(rkx) t
!
  m = mu + ml + 1
!
!  JOB = 0, Solve  a * x = b.
!
!  First solve l*y = b.
!
  if ( job == 0 ) then

     if ( ml > 0 ) then

        do k = 1, n-1
           lm = min ( ml, n-k )
           l = ipvt(k)
           t = b(l)
           if ( l /= k ) then
             b(l) = b(k)
             b(k) = t
           end if
           call saxpy ( lm, t, abd(m+1,k), 1, b(k+1), 1 )
        end do

     end if
!
!  Now solve u*x = y.
!
     do k = n, 1, -1
       b(k) = b(k) / abd(m,k)
       lm = min ( k, m ) - 1
       la = m - lm
       lb = k - lm
       t = -b(k)
       call saxpy ( lm, t, abd(la,k), 1, b(lb), 1 )
     end do
!
!  JOB nonzero, solve A' * X = B.
!
!  First solve U'*Y = B.
!
  else

     do k = 1, n
       lm = min ( k, m ) - 1
       t = dot_product ( abd(m-lm:m-1,k), b(k-lm:k-1) )
       b(k) = ( b(k) - t ) / abd(m,k)
     end do
!
!  Now solve L'*X = Y.
!
     if ( ml > 0 ) then

        do k = n-1, 1, -1
          lm = min ( ml, n-k )
          b(k) = b(k) + dot_product ( abd(m+1:m+lm,k), b(k+1:k+lm) )
          l = ipvt(k)
          if ( l /= k ) then
            call r_swap ( b(l), b(k) )
          end if
        end do

      end if

  end if

  return
end subroutine sgbsl

subroutine sgeco ( a, lda, n, ipvt, rcond, z )
!
!*******************************************************************************
!
!! SGECO factors a real(rkx) matrix and estimates its condition number.
!
!
!  Discussion:
!
!    If RCOND is not needed, SGEFA is slightly faster.
!
!    To solve A * X = B, follow SGECO by SGESL.
!
!    To compute inverse ( A ) * C, follow SGECO by SGESL.
!
!    To compute determinant ( A ), follow SGECO by SGEDI.
!
!    To compute inverse ( A ), follow SGECO by SGEDI.
!
!    For the system A * X = B, relative perturbations in A and B
!    of size EPSILON may cause relative perturbations in X of size
!    EPSILON/RCOND.
!
!    If RCOND is so small that the logical expression
!      1.0_rkx + RCOND == 1.0_rkx
!    is true, then A may be singular to working precision.  In particular,
!    RCOND is zero if exact singularity is detected or the estimate
!    underflows.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Author:
!
!    Cleve Moler,
!    University of New Mexico / Argonne National Lab.
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,N).  On input, a matrix to be factored.
!    On output, the LU factorization of the matrix.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix A.
!
!    Output, integer IPVT(N), the pivot indices.
!
!    Output, real(rkx) RCOND, an estimate of the reciprocal condition number of A.
!
!    Output, real(rkx) Z(N), a work vector whose contents are usually unimportant.
!    If A is close to a singular matrix, then Z is an approximate null vector
!    in the sense that
!      norm ( A * Z ) = RCOND * norm ( A ) * norm ( Z ).
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  real(rkx) anorm
  real(rkx) ek
  integer info
  integer ipvt(n)
  integer j
  integer k
  integer l
  real(rkx) rcond
  real(rkx) s
  real(rkx) sm
  real(rkx) t
  real(rkx) wk
  real(rkx) wkm
  real(rkx) ynorm
  real(rkx) z(n)
!
!  Compute the L1 norm of A.
!
  anorm = 0.0_rkx
  do j = 1, n
    anorm = max ( anorm, sum ( abs ( a(1:n,j) ) ) )
  end do
!
!  Compute the LU factorization.
!
  call sgefa ( a, lda, n, ipvt, info )
!
!  RCOND = 1 / ( norm(A) * (estimate of norm(inverse(A))) )
!
!  estimate of norm(inverse(A)) = norm(Z) / norm(Y)
!
!  where
!    A * Z = Y
!  and
!    A' * Y = E
!
!  The components of E are chosen to cause maximum local growth in the
!  elements of W, where U'*W = E.  The vectors are frequently rescaled
!  to avoid overflow.
!
!  Solve U' * W = E.
!
  ek = 1.0_rkx
  z(1:n) = 0.0_rkx

  do k = 1, n

    if ( z(k) /= 0.0_rkx ) then
      ek = sign ( ek, -z(k) )
    end if

    if ( abs ( ek - z(k) ) > abs ( a(k,k) ) ) then
      s = abs ( a(k,k) ) / abs ( ek - z(k) )
      z(1:n) = s * z(1:n)
      ek = s * ek
    end if

    wk = ek - z(k)
    wkm = -ek - z(k)
    s = abs ( wk )
    sm = abs ( wkm )

    if ( a(k,k) /= 0.0_rkx ) then
      wk = wk / a(k,k)
      wkm = wkm / a(k,k)
    else
      wk = 1.0_rkx
      wkm = 1.0_rkx
    end if

    if ( k+1 <= n ) then

      do j = k+1, n
        sm = sm + abs ( z(j) + wkm * a(k,j) )
        z(j) = z(j) + wk * a(k,j)
        s = s + abs ( z(j) )
      end do

      if ( s < sm ) then
        t = wkm - wk
        wk = wkm
        z(k+1:n) = z(k+1:n) + t * a(k,k+1:n)
      end if

    end if

    z(k) = wk

  end do

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )
!
!  Solve L' * Y = W
!
  do k = n, 1, -1

    z(k) = z(k) + dot_product ( a(k+1:n,k), z(k+1:n) )

    if ( abs ( z(k) ) > 1.0_rkx ) then
      z(1:n) = z(1:n) / abs ( z(k) )
    end if

    l = ipvt(k)

    t = z(l)
    z(l) = z(k)
    z(k) = t

  end do

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )

  ynorm = 1.0_rkx
!
!  Solve L * V = Y.
!
  do k = 1, n

    l = ipvt(k)

    t = z(l)
    z(l) = z(k)
    z(k) = t

    z(k+1:n) = z(k+1:n) + t * a(k+1:n,k)

    if ( abs ( z(k) ) > 1.0_rkx ) then
      ynorm = ynorm / abs ( z(k) )
      z(1:n) = z(1:n) / abs ( z(k) )
    end if

  end do

  s = sum ( abs ( z(1:n) ) )
  z(1:n) = z(1:n) / s
  ynorm = ynorm / s
!
!  Solve U * Z = V.
!
  do k = n, 1, -1

    if ( abs ( z(k) ) > abs ( a(k,k) ) ) then
      s = abs ( a(k,k) ) / abs ( z(k) )
      z(1:n) = s * z(1:n)
      ynorm = s * ynorm
    end if

    if ( a(k,k) /= 0.0_rkx ) then
      z(k) = z(k) / a(k,k)
    else
      z(k) = 1.0_rkx
    end if

    z(1:k-1) = z(1:k-1) - z(k) * a(1:k-1,k)

  end do
!
!  Normalize Z in the L1 norm.
!
  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = s * ynorm

  if ( anorm /= 0.0_rkx ) then
    rcond = ynorm / anorm
  else
    rcond = 0.0_rkx
  end if

  return
end subroutine sgeco

subroutine sgedi ( a, lda, n, ipvt, det, work, job )
!
!*******************************************************************************
!
!! SGEDI computes the determinant and inverse of a matrix factored by SGECO or SGEFA.
!
!
!  Discussion:
!
!    A division by zero will occur if the input factor contains
!    a zero on the diagonal and the inverse is requested.
!    It will not occur if the subroutines are called correctly
!    and if SGECO has set RCOND > 0.0_rkx or SGEFA has set INFO == 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,N), on input, the N by N factored matrix.
!    as output by SGECO or SGEFA.  On output, contains the inverse
!    matrix if requested.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix A.
!
!    Input, integer IPVT(N), the pivot vector from SGECO or SGEFA.
!
!    Workspace, real(rkx) WORK(N).
!
!    Output, real(rkx) DET(2), the determinant of original matrix if requested.
!    determinant = DET(1) * 10.0**DET(2)
!    with  1.0_rkx <= abs ( DET(1) ) < 10.0_rkx
!    or DET(1) == 0.0_rkx.
!
!    Input, integer JOB, specifies what is to be computed.
!    11, both determinant and inverse.
!    01, inverse only.
!    10, determinant only.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  real(rkx) det(2)
  integer i
  integer ipvt(n)
  integer j
  integer job
  integer k
  integer l
  real(rkx) t
  real(rkx), parameter :: ten = 10.0_rkx
  real(rkx) work(n)
!
!  Compute the determinant.
!
  if ( job / 10 /= 0 ) then

     det(1) = 1.0_rkx
     det(2) = 0.0_rkx

     do i = 1, n

        if ( ipvt(i) /= i ) then
          det(1) = - det(1)
        end if

        det(1) = a(i,i) * det(1)

        if ( det(1) == 0.0_rkx ) then
          exit
        end if

        do while ( abs ( det(1) ) < 1.0_rkx )
          det(1) = ten * det(1)
          det(2) = det(2) - 1.0_rkx
        end do

        do while ( abs ( det(1) ) >= ten )
          det(1) = det(1) / ten
          det(2) = det(2) + 1.0_rkx
        end do

    end do

  end if
!
!  Compute inverse(U).
!
  if ( mod ( job, 10 ) /= 0 ) then

     do k = 1, n

        a(k,k) = 1.0_rkx / a(k,k)
        t = - a(k,k)
        call sscal ( k-1, t, a(1,k), 1 )

        do j = k+1, n
           t = a(k,j)
           a(k,j) = 0.0_rkx
           call saxpy ( k, t, a(1,k), 1, a(1,j), 1 )
        end do

     end do
!
!  Form inverse(U) * inverse(L).
!
     do k = n-1, 1, -1

        do i = k+1, n
           work(i) = a(i,k)
           a(i,k) = 0.0_rkx
        end do

        do j = k+1, n
           t = work(j)
           call saxpy ( n, t, a(1,j), 1, a(1,k), 1 )
        end do

        l = ipvt(k)
        if ( l /= k ) then
          call sswap ( n, a(1,k), 1, a(1,l), 1 )
        end if

     end do

  end if

  return
end subroutine sgedi

subroutine sgefa ( a, lda, n, ipvt, info )
!
!*******************************************************************************
!
!! SGEFA factors a real(rkx) matrix.
!
!
!  Modified:
!
!    17 April 2002
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,N).
!    On intput, the matrix to be factored.
!    On output, an upper triangular matrix and the multipliers used to obtain
!    it.  The factorization can be written A=L*U, where L is a product of
!    permutation and unit lower triangular matrices, and U is upper triangular.
!
!    Input, integer LDA, the leading dimension of A.
!
!    Input, integer N, the order of the matrix A.
!
!    Output, integer IPVT(N), the pivot indices.
!
!    Output, integer INFO, singularity indicator.
!    0, normal value.
!    K, if U(K,K) == 0.  This is not an error condition for this subroutine,
!    but it does indicate that SGESL or SGEDI will divide by zero if called.
!    Use RCOND in SGECO for a reliable indication of singularity.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  integer info
  integer ipvt(n)
  integer j
  integer k
  integer l
  real(rkx) t
!
!  Gaussian elimination with partial pivoting.
!
  info = 0

  do k = 1, n - 1
!
!  Find L = pivot index.
!
    l = isamax ( n-k+1, a(k,k), 1 ) + k - 1
    ipvt(k) = l
!
!  Zero pivot implies this column already triangularized.
!
    if ( a(l,k) == 0.0_rkx ) then
      info = k
      cycle
    end if
!
!  Interchange if necessary.
!
    if ( l /= k ) then
      call r_swap ( a(l,k), a(k,k) )
    end if
!
!  Compute multipliers.
!
    a(k+1:n,k) = - a(k+1:n,k) / a(k,k)
!
!  Row elimination with column indexing.
!
    do j = k+1, n
      t = a(l,j)
      if ( l /= k ) then
        a(l,j) = a(k,j)
        a(k,j) = t
      end if
      call saxpy ( n-k, t, a(k+1,k), 1, a(k+1,j), 1 )
    end do

  end do

  ipvt(n) = n

  if ( a(n,n) == 0.0_rkx ) then
    info = n
  end if

  return
end subroutine sgefa

subroutine sgesl ( a, lda, n, ipvt, b, job )
!
!*******************************************************************************
!
!! SGESL solves a real(rkx) general linear system A * X = B.
!
!
!  Discussion:
!
!    SGESL can solve either of the systems A * X = B or A' * X = B.
!
!    The system matrix must have been factored by SGECO or SGEFA.
!
!    A division by zero will occur if the input factor contains a
!    zero on the diagonal.  Technically this indicates singularity
!    but it is often caused by improper arguments or improper
!    setting of LDA.  It will not occur if the subroutines are
!    called correctly and if SGECO has set RCOND > 0.0_rkx
!    or SGEFA has set INFO == 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Modified:
!
!    17 April 2002
!
!  Parameters:
!
!    Input, real(rkx) A(LDA,N), the output from SGECO or SGEFA.
!
!    Input, integer LDA, the leading dimension of A.
!
!    Input, integer N, the order of the matrix A.
!
!    Input, integer IPVT(N), the pivot vector from SGECO or SGEFA.
!
!    Input/output, real(rkx) B(N).
!    On input, the right hand side vector.
!    On output, the solution vector.
!
!    Input, integer JOB.
!    0, solve A * X = B;
!    nonzero, solve A' * X = B.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  real(rkx) b(n)
  integer ipvt(n)
  integer job
  integer k
  integer l
  real(rkx) t
!
!  Solve A * X = B.
!
  if ( job == 0 ) then

    do k = 1, n-1

      l = ipvt(k)
      t = b(l)

      if ( l /= k ) then
        b(l) = b(k)
        b(k) = t
      end if

      call saxpy ( n-k, t, a(k+1,k), 1, b(k+1), 1 )

    end do

    do k = n, 1, -1
      b(k) = b(k) / a(k,k)
      t = -b(k)
      call saxpy ( k-1, t, a(1,k), 1, b(1), 1 )
    end do

  else
!
!  Solve A' * X = B.
!
    do k = 1, n
      t = dot_product ( a(1:k-1,k), b(1:k-1) )
      b(k) = ( b(k) - t ) / a(k,k)
    end do

    do k = n-1, 1, -1

      b(k) = b(k) + dot_product ( a(k+1:n,k), b(k+1:n) )
      l = ipvt(k)

      if ( l /= k ) then
        call r_swap ( b(l), b(k) )
      end if

    end do

  end if

  return
end subroutine sgesl

subroutine sgtsl ( n, c, d, e, b, info )
!
!*******************************************************************************
!
!! SGTSL solves a general tridiagonal linear system.
!
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Modified:
!
!    31 October 2001
!
!  Parameters:
!
!    Input, integer N, the order of the tridiagonal matrix.
!
!    Input/output, real(rkx) C(N), contains the subdiagonal of the tridiagonal
!    matrix in entries C(2:N).  On output, C is destroyed.
!
!    Input/output, real(rkx) D(N).  On input, the diagonal of the matrix.
!    On output, D is destroyed.
!
!    Input/output, real(rkx) E(N), contains the superdiagonal of the tridiagonal
!    matrix in entries E(1:N-1).  On output E is destroyed.
!
!    Input/output, real(rkx) B(N).  On input, the right hand side.  On output,
!    the solution.
!
!    Output, integer INFO, error flag.
!    0, normal value.
!    K, the K-th element of the diagonal becomes exactly zero.  The
!       subroutine returns if this error condition is detected.
!
  implicit none
!
  integer n
!
  real(rkx) b(n)
  real(rkx) c(n)
  real(rkx) d(n)
  real(rkx) e(n)
  integer info
  integer k
  real(rkx) t
!
  info = 0
  c(1) = d(1)

  if ( n >= 2 ) then

    d(1) = e(1)
    e(1) = 0.0_rkx
    e(n) = 0.0_rkx

    do k = 1, n - 1
!
!  Find the larger of the two rows, and interchange if necessary.
!
      if ( abs ( c(k+1) ) >= abs ( c(k) ) ) then
        call r_swap ( c(k), c(k+1) )
        call r_swap ( d(k), d(k+1) )
        call r_swap ( e(k), e(k+1) )
        call r_swap ( b(k), b(k+1) )
      end if
!
!  Fail if no nonzero pivot could be found.
!
      if ( c(k) == 0.0_rkx ) then
        info = k
        return
      end if
!
!  Zero elements.
!
      t = -c(k+1) / c(k)
      c(k+1) = d(k+1) + t * d(k)
      d(k+1) = e(k+1) + t * e(k)
      e(k+1) = 0.0_rkx
      b(k+1) = b(k+1) + t * b(k)

    end do

  end if

  if ( c(n) == 0.0_rkx ) then
    info = n
    return
  end if
!
!  Back solve.
!
  b(n) = b(n) / c(n)

  if ( n > 1 ) then

    b(n-1) = ( b(n-1) - d(n-1) * b(n) ) / c(n-1)

    do k = n-2, 1, -1
      b(k) = ( b(k) - d(k) * b(k+1) - e(k) * b(k+2) ) / c(k)
    end do

  end if

  return
end subroutine sgtsl

function snrm2 ( n, x, incx )
!
!*******************************************************************************
!
!! SNRM2 computes the Euclidean norm of a vector.
!
!
!  Discussion:
!
!    The original SNRM2 algorithm is accurate but written in an
!    obscure and obsolete format.  This version goes for clarity.
!
!  Modified:
!
!    01 June 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real(rkx) X(*), the vector whose norm is to be computed.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real(rkx) SNRM2, the Euclidean norm of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer n
  real(rkx) snrm2
  real(rkx) stemp
  real(rkx) x(*)
  real(rkx) xmax
!
  if ( n <= 0 ) then

    snrm2 = 0.0_rkx

  else

    xmax = samax ( n, x, incx )

    if ( xmax == 0.0_rkx ) then

      snrm2 = 0.0_rkx

    else

      if ( incx >= 0 ) then
        ix = 1
      else
        ix = ( - n + 1 ) * incx + 1
      end if

      stemp = 0.0_rkx
      do i = 1, n
        stemp = stemp + ( x(ix) / xmax )**2
        ix = ix + incx
      end do

      snrm2 = xmax * sqrt ( stemp )

    end if

  end if

  return
end function snrm2

subroutine spbco ( abd, lda, n, m, rcond, z, info )
!
!*******************************************************************************
!
!! SPBCO factors a real(rkx) symmetric positive definite banded matrix.
!
!
!  Discussion:
!
!    SPBCO also estimates the condition of the matrix.
!
!    If RCOND is not needed, SPBFA is slightly faster.
!
!    To solve A*X = B, follow SPBCO by SPBSL.
!
!    To compute inverse(A)*C, follow SPBCO by SPBSL.
!
!    To compute determinant(A), follow SPBCO by SPBDI.
!
!  Band storage:
!
!    If A is a symmetric positive definite band matrix, the following
!    program segment will set up the input.
!
!      m = (band width above diagonal)
!      do j = 1, n
!        i1 = max (1, j-m)
!        do i = i1, j
!          k = i-j+m+1
!          abd(k,j) = a(i,j)
!        end do
!      end do
!
!    This uses M + 1 rows of A, except for the M by M upper left triangle,
!    which is ignored.
!
!    For example, if the original matrix is
!
!      11 12 13  0  0  0
!      12 22 23 24  0  0
!      13 23 33 34 35  0
!       0 24 34 44 45 46
!       0  0 35 45 55 56
!       0  0  0 46 56 66
!
!    then N = 6, M = 2  and ABD should contain
!
!       *  * 13 24 35 46
!       * 12 23 34 45 56
!      11 22 33 44 55 66
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) ABD(LDA,N).  On input, the matrix to be factored.
!    The columns of the upper triangle are stored in the columns of ABD and the
!    diagonals of the upper triangle are stored in the rows of ABD.
!    On output, an upper triangular matrix R, stored in band form, so that
!    A = R'*R.  If INFO /= 0, the factorization is not complete.
!
!    Input, integer LDA, the leading dimension of the array ABD.
!    LDA must be >= M+1.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer M, the number of diagonals above the main diagonal.
!
!    Output, real(rkx) RCOND, an estimate of the reciprocal condition of A.
!    For the system A*X = B, relative perturbations in A and B of size
!    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
!    If RCOND is so small that the logical expression
!      1.0_rkx + RCOND == 1.0_rkx
!    is true, then A may be singular to working precision.  In particular,
!    RCOND is zero if exact singularity is detected or the estimate underflows.
!
!    Output, real(rkx) Z(N), a work vector whose contents are usually unimportant.
!    If A is singular to working precision, then Z is an approximate null
!    vector in the sense that
!      norm(A*Z) = RCOND * norm(A) * norm(Z).
!    If INFO /= 0, Z is unchanged.
!
!    Output, integer INFO, error flag.
!    0, for normal return.
!    K, signals an error condition.  The leading minor of order K is not
!    positive definite.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) abd(lda,n)
  real(rkx) anorm
  real(rkx) ek
  integer i
  integer info
  integer j
  integer j2
  integer k
  integer l
  integer la
  integer lb
  integer lm
  integer m
  integer mu
  real(rkx) rcond
  real(rkx) s
  real(rkx) sm
  real(rkx) t
  real(rkx) wk
  real(rkx) wkm
  real(rkx) ynorm
  real(rkx) z(n)
!
!  Find the norm of A.
!
  do j = 1, n

    l = min ( j, m+1 )
    mu = max ( m+2-j, 1 )
    z(j) = sum ( abs ( abd(mu:mu+l-1,j) ) )
    k = j - l
    do i = mu, m
      k = k + 1
      z(k) = z(k) + abs ( abd(i,j) )
    end do

  end do

  anorm = maxval ( z(1:n) )
!
!  Factor.
!
  call spbfa ( abd, lda, n, m, info )

  if ( info /= 0 ) then
    return
  end if
!
!  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
!
!  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
!
!  The components of E are chosen to cause maximum local
!  growth in the elements of W where R'*W = E.
!
!  The vectors are frequently rescaled to avoid overflow.
!
!  Solve R'*W = E.
!
  ek = 1.0_rkx
  z(1:n) = 0.0_rkx

  do k = 1, n

    if ( z(k) /= 0.0_rkx ) then
      ek = sign ( ek, -z(k) )
    end if

    if ( abs ( ek - z(k) ) > abd(m+1,k) ) then
      s = abd(m+1,k) / abs ( ek - z(k) )
      z(1:n) = s * z(1:n)
      ek = s * ek
    end if

    wk = ek - z(k)
    wkm = -ek - z(k)
    s = abs ( wk )
    sm = abs ( wkm )
    wk = wk / abd(m+1,k)
    wkm = wkm / abd(m+1,k)
    j2 = min ( k+m, n )
    i = m + 1

    if ( k+1 <= j2 ) then

      do j = k+1, j2
        i = i - 1
        sm = sm + abs ( z(j) + wkm * abd(i,j) )
        z(j) = z(j) + wk * abd(i,j)
        s = s + abs ( z(j) )
      end do

      if ( s < sm ) then

        t = wkm - wk
        wk = wkm
        i = m + 1

        do j = k+1, j2
          i = i - 1
          z(j) = z(j) + t * abd(i,j)
        end do

      end if

    end if

    z(k) = wk

  end do

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )
!
!  Solve R*Y = W.
!
  do k = n, 1, -1

    if ( abs ( z(k) ) > abd(m+1,k) ) then
      s = abd(m+1,k) / abs ( z(k) )
      z(1:n) = s * z(1:n)
    end if

    z(k) = z(k) / abd(m+1,k)
    lm = min ( k-1, m )
    la = m + 1 - lm
    lb = k - lm
    t = -z(k)
    call saxpy ( lm, t, abd(la,k), 1, z(lb), 1 )

  end do

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )

  ynorm = 1.0_rkx
!
!  Solve R'*V = Y.
!
  do k = 1, n

    lm = min ( k-1, m )
    la = m + 1 - lm
    lb = k - lm

    z(k) = z(k) - dot_product ( abd(la:la+lm-1,k), z(lb:lb+lm-1) )

    if ( abs ( z(k) ) > abd(m+1,k) ) then
      s = abd(m+1,k) / abs ( z(k) )
      z(1:n) = s * z(1:n)
      ynorm = s * ynorm
    end if

    z(k) = z(k) / abd(m+1,k)

  end do

  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = s * ynorm
!
!  Solve R*Z = W.
!
  do k = n, 1, -1

    if ( abs ( z(k) ) > abd(m+1,k) ) then
      s = abd(m+1,k) / abs ( z(k) )
      z(1:n) = s * z(1:n)
      ynorm = s * ynorm
    end if

    z(k) = z(k) / abd(m+1,k)
    lm = min ( k-1, m )
    la = m + 1 - lm
    lb = k - lm
    t = -z(k)
    call saxpy ( lm, t, abd(la,k), 1, z(lb), 1 )

  end do
!
!  Make ZNORM = 1.0.
!
  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = s * ynorm

  if ( anorm /= 0.0_rkx ) then
    rcond = ynorm / anorm
  else
    rcond = 0.0_rkx
  end if

  return
end subroutine spbco

subroutine spbdi ( abd, lda, n, m, det )
!
!*******************************************************************************
!
!! SPBDI computes the determinant of a matrix factored by SPBCO or SPBFA.
!
!
!  Discussion:
!
!    If the inverse is needed, use SPBSL N times.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) ABD(LDA,N), the output from SPBCO or SPBFA.
!
!    Input, integer LDA, the leading dimension of the array ABD.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer M, the number of diagonals above the main diagonal.
!
!    Output, real(rkx) DET(2), the determinant of the original matrix in the form
!      determinant = DET(1) * 10.0**DET(2)
!    with 1.0_rkx <= DET(1) < 10.0_rkx or DET(1) == 0.0_rkx.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) abd(lda,n)
  real(rkx) det(2)
  integer i
  integer m
  real(rkx) s
!
!  Compute the determinant.
!
  det(1) = 1.0_rkx
  det(2) = 0.0_rkx
  s = 10.0_rkx

  do i = 1, n

    det(1) = abd(m+1,i)**2 * det(1)

    if ( det(1) == 0.0_rkx ) then
      return
    end if

    do while ( det(1) < 1.0_rkx )
      det(1) = s * det(1)
      det(2) = det(2) - 1.0_rkx
    end do

    do while ( det(1) >= s )
      det(1) = det(1) / s
      det(2) = det(2) + 1.0_rkx
    end do

  end do

  return
end subroutine spbdi

subroutine spbfa ( abd, lda, n, m, info )
!
!*******************************************************************************
!
!! SPBFA factors a real(rkx) symmetric positive definite matrix stored in band form.
!
!
!  Discussion:
!
!    SPBFA is usually called by SPBCO, but it can be called
!    directly with a saving in time if RCOND is not needed.
!
!  Band storage:
!
!    If A is a symmetric positive definite band matrix,
!    the following program segment will set up the input.
!
!      m = (band width above diagonal)
!      do j = 1, n
!        i1 = max (1, j-m)
!        do i = i1, j
!          k = i-j+m+1
!          abd(k,j) = a(i,j)
!        end do
!      end do
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) ABD(LDA,N).  On input, the matrix to be factored.
!    The columns of the upper triangle are stored in the columns of ABD
!    and the diagonals of the upper triangle are stored in the
!    rows of ABD.  On output, an upper triangular matrix R, stored in band
!    form, so that A = R' * R.
!
!    Input, integer LDA, the leading dimension of the array ABD.
!    LDA must be >= M+1.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer M, the number of diagonals above the main diagonal.
!
!    Output, integer INFO, error indicator.
!    0, for normal return.
!    K, if the leading minor of order K is not positive definite.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) abd(lda,n)
  integer ik
  integer info
  integer j
  integer jk
  integer k
  integer m
  integer mu
  real(rkx) s
  real(rkx) t
!
  do j = 1, n

    s = 0.0_rkx
    ik = m + 1
    jk = max (j-m,1)
    mu = max (m+2-j,1)

    do k = mu, m
      t = abd(k,j) - sdot ( k-mu, abd(ik,jk), 1, abd(mu,j), 1 )
      t = t / abd(m+1,jk)
      abd(k,j) = t
      s = s + t * t
      ik = ik - 1
      jk = jk + 1
    end do

    s = abd(m+1,j) - s

    if ( s <= 0.0_rkx ) then
      info = j
      return
    end if

    abd(m+1,j) = sqrt ( s )

  end do

  info = 0

  return
end subroutine spbfa

subroutine spbsl ( abd, lda, n, m, b )
!
!*******************************************************************************
!
!! SPBSL solves a real(rkx) symmetric positive definite band factored by SPBCO or SPBFA.
!
!
!  Discussion:
!
!    To compute inverse(A) * C  where C is a matrix with P columns:
!
!      call spbco ( abd, lda, n, rcond, z, info )
!
!      if ( rcond is too small .or. info /= 0) go to ...
!
!      do j = 1, p
!        call spbsl ( abd, lda, n, c(1,j) )
!      end do
!
!    A division by zero will occur if the input factor contains
!    a zero on the diagonal.  Technically this indicates
!    singularity but it is usually caused by improper subroutine
!    arguments.  It will not occur if the subroutines are called
!    correctly and INFO == 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) ABD(LDA,N), the output from SPBCO or SPBFA.
!
!    Input, integer LDA, the leading dimension of the array ABD.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer M, the number of diagonals above the main diagonal.
!
!    Input/output, real(rkx) B(N).  On input, the right hand side.
!    On output, the solution.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) abd(lda,n)
  real(rkx) b(n)
  integer k
  integer la
  integer lb
  integer lm
  integer m
  real(rkx) t
!
!  Solve R'*Y = B.
!
  do k = 1, n
     lm = min ( k-1, m )
     t = dot_product ( abd(m+1-lm:m,k), b(k-lm:k-1) )
     b(k) = ( b(k) - t ) / abd(m+1,k)
  end do
!
!  Solve R*X = Y.
!
  do k = n, 1, -1
     lm = min ( k-1, m )
     la = m + 1 - lm
     lb = k - lm
     b(k) = b(k) / abd(m+1,k)
     t = -b(k)
     call saxpy ( lm, t, abd(la,k), 1, b(lb), 1 )
  end do

  return
end subroutine spbsl

subroutine spoco ( a, lda, n, rcond, z, info )
!
!*******************************************************************************
!
!! SPOCO factors a real(rkx) symmetric positive definite matrix and estimates its condition.
!
!
!  Discussion:
!
!     If RCOND is not needed, SPOFA is slightly faster.
!
!     To solve A*X = B, follow SPOCO by SPOSL.
!
!     To compute inverse(A)*C, follow SPOCO by SPOSL.
!
!     To compute determinant(A), follow SPOCO by SPODI.
!
!     To compute inverse(A), follow SPOCO by SPODI.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,N).  On input, the symmetric matrix to
!    be factored.  Only the diagonal and upper triangle are used.
!    On output, an upper triangular matrix R so that A = R'*R where R'
!    is the transpose.  The strict lower triangle is unaltered.
!    If INFO /= 0, the factorization is not complete.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix.
!
!    Output, real(rkx) RCOND, an estimate of the reciprocal condition of A.
!    For the system A*X = B, relative perturbations in A and B of size
!    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
!    If RCOND is so small that the logical expression
!      1.0_rkx + RCOND == 1.0_rkx
!    is true, then A may be singular to working precision.  In particular,
!    RCOND is zero if exact singularity is detected or the estimate underflows.
!
!    Output, real(rkx) Z(N), a work vector whose contents are usually unimportant.
!    If A is close to a singular matrix, then Z is an approximate null vector
!    in the sense that
!      norm(A*Z) = RCOND * norm(A) * norm(Z).
!    If INFO /= 0, Z is unchanged.
!
!    Output, integer INFO, error flag.
!    0, for normal return.
!    K, signals an error condition.  The leading minor of order K is not
!    positive definite.
!
  implicit none
!
  integer lda
  integer n

  real(rkx) a(lda,n)
  real(rkx) anorm
  real(rkx) ek
  integer i
  integer info
  integer j
  integer k
  real(rkx) rcond
  real(rkx) s
  real(rkx) sm
  real(rkx) t
  real(rkx) wk
  real(rkx) wkm
  real(rkx) ynorm
  real(rkx) z(n)
!
!  Find norm of A using only upper half.
!
  do j = 1, n
    z(j) = sum ( abs ( a(1:j,j) ) )
    do i = 1, j-1
      z(i) = z(i) + abs ( a(i,j) )
    end do
  end do

  anorm = maxval ( z(1:n) )
!
!  Factor.
!
  call spofa ( a, lda, n, info )

  if ( info /= 0 ) then
    return
  end if
!
!  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
!
!  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
!
!  The components of E are chosen to cause maximum local
!  growth in the elements of W where R'*W = E.
!
!  The vectors are frequently rescaled to avoid overflow.
!
!  Solve R'*W = E.
!
     ek = 1.0_rkx
     z(1:n) = 0.0_rkx

     do k = 1, n

        if ( z(k) /= 0.0_rkx ) then
          ek = sign ( ek, -z(k) )
        end if

        if ( abs ( ek - z(k) ) > a(k,k) ) then
          s = a(k,k) / abs ( ek - z(k) )
          z(1:n) = s * z(1:n)
          ek = s * ek
        end if

        wk = ek - z(k)
        wkm = -ek - z(k)
        s = abs ( wk )
        sm = abs ( wkm )
        wk = wk / a(k,k)
        wkm = wkm / a(k,k)

        if ( k+1 <= n ) then

           do j = k+1, n
              sm = sm + abs ( z(j) + wkm * a(k,j) )
              z(j) = z(j) + wk * a(k,j)
              s = s + abs ( z(j) )
           end do

           if ( s < sm ) then

              t = wkm - wk
              wk = wkm
              do j = k+1, n
                z(j) = z(j) + t * a(k,j)
              end do

           end if

        end if

        z(k) = wk

     end do

     z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )
!
!  Solve R*Y = W.
!
     do k = n, 1, -1

        if ( abs ( z(k) ) > a(k,k) ) then
           s = a(k,k) / abs ( z(k) )
           z(1:n) = s * z(1:n)
        end if

        z(k) = z(k)/a(k,k)
        t = -z(k)
        call saxpy ( k-1, t, a(1,k), 1, z(1), 1 )

     end do

     z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )
     ynorm = 1.0_rkx
!
!  Solve R'*V = Y.
!
     do k = 1, n

        z(k) = z(k) - dot_product ( a(1:k-1,k), z(1:k-1) )

        if ( abs ( z(k) ) > a(k,k) ) then
           s = a(k,k) / abs ( z(k) )
           z(1:n) = s * z(1:n)
           ynorm = s * ynorm
        end if

        z(k) = z(k) / a(k,k)

     end do

     s = 1.0_rkx / sum ( abs ( z(1:n) ) )
     z(1:n) = s * z(1:n)
     ynorm = s * ynorm
!
!  Solve R*Z = V.
!
     do k = n, 1, -1

        if ( abs ( z(k) ) > a(k,k) ) then
           s = a(k,k) / abs ( z(k) )
           z(1:n) = s * z(1:n)
           ynorm = s * ynorm
        end if

        z(k) = z(k)/a(k,k)
        t = -z(k)
        call saxpy ( k-1, t, a(1,k), 1, z(1), 1 )

     end do
!
!  Make ZNORM = 1.0.
!
     s = 1.0_rkx / sum ( abs ( z(1:n) ) )
     z(1:n) = s * z(1:n)
     ynorm = s * ynorm

     if ( anorm /= 0.0_rkx ) then
       rcond = ynorm / anorm
     else
       rcond = 0.0_rkx
     end if

  return
end subroutine spoco

subroutine spodi ( a, lda, n, det, job )
!
!*******************************************************************************
!
!! SPODI computes the determinant and inverse of a certain matrix.
!
!
!  Discussion:
!
!    The matrix is real(rkx) symmetric positive definite.
!    SPODI uses the factors computed by SPOCO, SPOFA or SQRDC.
!
!    A division by zero will occur if the input factor contains
!    a zero on the diagonal and the inverse is requested.
!    It will not occur if the subroutines are called correctly
!    and if SPOCO or SPOFA has set INFO == 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,N).  On input, the output A from SPOCO
!    or SPOFA, or the output X from SQRDC.  On output, if SPOCO or
!    SPOFA was used to factor A then SPODI produces the upper half of
!    inverse(A).  If SQRDC was used to decompose X then SPODI produces
!    the upper half of inverse(X'*X) where X' is the transpose.
!    Elements of A below the diagonal are unchanged.  If the units digit
!    of JOB is zero, A is unchanged.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix A.
!
!    Input, integer JOB, specifies the task.
!    11, both determinant and inverse.
!    01, inverse only.
!    10, determinant only.
!
!    Output, real(rkx) DET(2), the determinant of A or of X'*X if requested.
!      determinant = DET(1) * 10.0**DET(2)
!    with 1.0_rkx <= DET(1) < 10.0_rkx or DET(1) == 0.0_rkx.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  real(rkx) det(2)
  integer i
  integer j
  integer job
  integer k
  real(rkx) s
  real(rkx) t
!
!  Compute the determinant.
!
  if ( job / 10 /= 0 ) then

     det(1) = 1.0_rkx
     det(2) = 0.0_rkx
     s = 10.0_rkx

     do i = 1, n

        det(1) = a(i,i)**2 * det(1)

        if ( det(1) == 0.0_rkx ) then
          exit
        end if

        do while ( det(1) < 1.0_rkx )
          det(1) = s * det(1)
          det(2) = det(2) - 1.0_rkx
        end do

        do while ( det(1) >= s )
          det(1) = det(1) / s
          det(2) = det(2) + 1.0_rkx
        end do

    end do

  end if
!
!  Compute inverse(R).
!
  if ( mod ( job, 10 ) /= 0 ) then

     do k = 1, n

        a(k,k) = 1.0_rkx / a(k,k)
        t = -a(k,k)
        call sscal ( k-1, t, a(1,k), 1 )

        do j = k+1, n
          t = a(k,j)
          a(k,j) = 0.0_rkx
          call saxpy ( k, t, a(1,k), 1, a(1,j), 1 )
        end do

     end do
!
!  Form inverse(R) * transpose(inverse(R)).
!
     do j = 1, n
        do k = 1, j-1
          t = a(k,j)
          call saxpy ( k, t, a(1,j), 1, a(1,k), 1 )
        end do
        t = a(j,j)
        call sscal ( j, t, a(1,j), 1 )
     end do

  end if

  return
end subroutine spodi

subroutine spofa ( a, lda, n, info )
!
!*******************************************************************************
!
!! SPOFA factors a real(rkx) symmetric positive definite matrix.
!
!
!  Discussion:
!
!    SPOFA is usually called by SPOCO, but it can be called
!    directly with a saving in time if RCOND is not needed.
!
!  Modified:
!
!    17 April 2002
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,N).  On input, the symmetric matrix to be
!    factored.  Only the diagonal and upper triangle are used.
!    On output, an upper triangular matrix R so that A = R'*R
!    where R' is the transpose.  The strict lower triangle is unaltered.
!    If INFO /= 0, the factorization is not complete.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix.
!
!    Output, integer INFO, error flag.
!    0, for normal return.
!    K, signals an error condition.  The leading minor of order K is not
!    positive definite.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  integer info
  integer j
  integer k
  real(rkx) s
  real(rkx) t
!
  do j = 1, n

    s = 0.0_rkx

    do k = 1, j-1
      t = a(k,j) - dot_product ( a(1:k-1,k), a(1:k-1,j) )
      t = t / a(k,k)
      a(k,j) = t
      s = s + t * t
    end do

    s = a(j,j) - s

    if ( s <= 0.0_rkx ) then
      info = j
      return
    end if

    a(j,j) = sqrt ( s )

  end do

  info = 0

  return
end subroutine spofa

subroutine sposl ( a, lda, n, b )
!
!*******************************************************************************
!
!! SPOSL solves a linear system factored by SPOCO or SPOFA.
!
!
!  Discussion:
!
!    To compute inverse(A) * C where C is a matrix with P columns:
!
!      call spoco ( a, lda, n, rcond, z, info )
!
!      if ( rcond is not too small .and. info == 0 ) then
!        do j = 1, p
!          call sposl ( a, lda, n, c(1,j) )
!        end do
!      end if
!
!    A division by zero will occur if the input factor contains
!    a zero on the diagonal.  Technically this indicates
!    singularity but it is usually caused by improper subroutine
!    arguments.  It will not occur if the subroutines are called
!    correctly and INFO == 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) A(LDA,N), the output from SPOCO or SPOFA.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix.
!
!    Input/output, real(rkx) B(N).  On input, the right hand side.
!    On output, the solution.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  real(rkx) b(n)
  integer k
  real(rkx) t
!
!  Solve R'*Y = B.
!
  do k = 1, n
    t = dot_product ( a(1:k-1,k), b(1:k-1) )
    b(k) = ( b(k) - t ) / a(k,k)
  end do
!
!  Solve R*X = Y.
!
  do k = n, 1, -1
    b(k) = b(k)/a(k,k)
    t = -b(k)
    call saxpy ( k-1, t, a(1,k), 1, b(1), 1 )
  end do

  return
end subroutine sposl

subroutine sppco ( ap, n, rcond, z, info )
!
!*******************************************************************************
!
!! SPPCO factors a real(rkx) symmetric positive definite matrix in packed form.
!
!
!  Discussion:
!
!    SPPCO also estimates the condition of the matrix.
!
!    If RCOND is not needed, SPPFA is slightly faster.
!
!    To solve A*X = B, follow SPPCO by SPPSL.
!
!    To compute inverse(A)*C, follow SPPCO by SPPSL.
!
!    To compute determinant(A), follow SPPCO by SPPDI.
!
!    To compute inverse(A), follow SPPCO by SPPDI.
!
!  Packed storage:
!
!    The following program segment will pack the upper triangle of
!    a symmetric matrix.
!
!      k = 0
!      do j = 1, n
!        do i = 1, j
!          k = k + 1
!          ap(k) = a(i,j)
!        end do
!      end do
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) AP(N*(N+1)/2).  On input, the packed form of a
!    symmetric matrix A.  The columns of the upper triangle are stored
!    sequentially in a one-dimensional array.  On output, an upper
!    triangular matrix R, stored in packed form, so that A = R'*R.
!    If INFO /= 0, the factorization is not complete.
!
!    Input, integer N, the order of the matrix.
!
!    Output, real(rkx) RCOND, an estimate of the reciprocal condition of A.
!    For the system A*X = B, relative perturbations in A and B of size
!    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
!    If RCOND is so small that the logical expression
!      1.0_rkx + RCOND == 1.0_rkx
!    is true, then A may be singular to working precision.  In particular,
!    RCOND is zero if exact singularity is detected or the estimate underflows.
!
!    Output, real(rkx) Z(N), a work vector whose contents are usually unimportant.
!    If A is singular to working precision, then Z is an approximate null
!    vector in the sense that
!      norm(A*Z) = RCOND * norm(A) * norm(Z).
!    If INFO /= 0, Z is unchanged.
!
!    Output, integer INFO, error flag.
!    0, for normal return.
!    K, signals an error condition.  The leading minor of order K is
!    not positive definite.
!
  implicit none
!
  integer n
!
  real(rkx) anorm
  real(rkx) ap((n*(n+1))/2)
  real(rkx) ek
  integer i
  integer ij
  integer info
  integer j
  integer j1
  integer k
  integer kj
  integer kk
  real(rkx) rcond
  real(rkx) s
  real(rkx) sm
  real(rkx) t
  real(rkx) wk
  real(rkx) wkm
  real(rkx) ynorm
  real(rkx) z(n)
!
!  Find the norm of A.
!
  j1 = 1
  do j = 1, n
    z(j) = sum ( abs ( ap(j1:j1+j-1) ) )
    ij = j1
    j1 = j1 + j
    do i = 1, j-1
      z(i) = z(i) + abs ( ap(ij) )
      ij = ij + 1
    end do
  end do

  anorm = maxval ( z(1:n) )
!
!  Factor.
!
  call sppfa ( ap, n, info )

  if ( info /= 0 ) then
    return
  end if
!
!  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
!
!  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
!
!  The components of E are chosen to cause maximum local
!  growth in the elements of W where R'*W = E.
!
!  The vectors are frequently rescaled to avoid overflow.
!
!  Solve R'*W = E.
!
     ek = 1.0_rkx
     z(1:n) = 0.0_rkx

     kk = 0

     do k = 1, n

        kk = kk + k

        if ( z(k) /= 0.0_rkx ) then
          ek = sign ( ek, -z(k) )
        end if

        if ( abs ( ek - z(k) ) > ap(kk) ) then
           s = ap(kk) / abs ( ek - z(k) )
           z(1:n) = s * z(1:n)
           ek = s * ek
        end if

        wk = ek - z(k)
        wkm = -ek - z(k)
        s = abs ( wk )
        sm = abs ( wkm )
        wk = wk / ap(kk)
        wkm = wkm / ap(kk)
        kj = kk + k

        if ( k+1 <= n ) then

           do j = k+1, n
              sm = sm + abs ( z(j) + wkm * ap(kj) )
              z(j) = z(j) + wk * ap(kj)
              s = s + abs ( z(j) )
              kj = kj + j
           end do

           if ( s < sm ) then

              t = wkm - wk
              wk = wkm
              kj = kk + k

              do j = k+1, n
                 z(j) = z(j) + t * ap(kj)
                 kj = kj + j
              end do

           end if

        end if

        z(k) = wk

     end do

     z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )
!
!  Solve R*Y = W.
!
     do k = n, 1, -1

        if ( abs ( z(k) ) > ap(kk) ) then
          s = ap(kk) / abs ( z(k) )
          z(1:n) = s * z(1:n)
        end if

        z(k) = z(k) / ap(kk)
        kk = kk - k
        t = -z(k)
        call saxpy ( k-1, t, ap(kk+1), 1, z(1), 1 )

     end do

     z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )

     ynorm = 1.0_rkx
!
!  Solve R'*V = Y.
!
     do k = 1, n

        z(k) = z(k) - sdot ( k-1, ap(kk+1), 1, z(1), 1 )
        kk = kk + k

        if ( abs ( z(k) ) > ap(kk) ) then
          s = ap(kk) / abs ( z(k) )
          z(1:n) = s * z(1:n)
          ynorm = s * ynorm
        end if

        z(k) = z(k) / ap(kk)

     end do

     s = 1.0_rkx / sum ( abs ( z(1:n) ) )
     z(1:n) = s * z(1:n)
     ynorm = s * ynorm
!
!  Solve R*Z = V.
!
     do k = n, 1, -1

        if ( abs ( z(k) ) > ap(kk) ) then
           s = ap(kk) / abs ( z(k) )
           z(1:n) = s * z(1:n)
           ynorm = s * ynorm
        end if

        z(k) = z(k) / ap(kk)
        kk = kk - k
        t = -z(k)
        call saxpy ( k-1, t, ap(kk+1), 1, z(1), 1 )

     end do
!
!  Make ZNORM = 1.0.
!
     s = 1.0_rkx / sum ( abs ( z(1:n) ) )
     z(1:n) = s * z(1:n)
     ynorm = s * ynorm

     if ( anorm /= 0.0_rkx ) then
       rcond = ynorm / anorm
     else
       rcond = 0.0_rkx
     end if

  return
end subroutine sppco

subroutine sppdi ( ap, n, det, job )
!
!*******************************************************************************
!
!! SPPDI computes the determinant and inverse of a matrix factored by SPPCO or SPPFA.
!
!
!  Discussion:
!
!    A division by zero will occur if the input factor contains
!    a zero on the diagonal and the inverse is requested.
!    It will not occur if the subroutines are called correctly
!    and if SPOCO or SPOFA has set INFO == 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) AP(N*(N+1)/2).  On input, the output from
!    SPPCO or SPPFA.  On output, the upper triangular half of the
!    inverse, if requested.
!
!    Input, integer N, the order of the matrix.
!
!    Output, real(rkx) DET(2), the determinant of the original matrix if requested.
!      determinant = DET(1) * 10.0**DET(2)
!    with  1.0_rkx <= DET(1) < 10.0_rkx or DET(1) == 0.0_rkx.
!
!    Input, integer JOB, job request.
!    11, both determinant and inverse.
!    01, inverse only.
!    10, determinant only.
!
  implicit none
!
  integer n
!
  real(rkx) ap((n*(n+1))/2)
  real(rkx) det(2)
  integer i
  integer ii
  integer j
  integer j1
  integer jj
  integer job
  integer k
  integer k1
  integer kj
  integer kk
  real(rkx) s
  real(rkx) t
!
!  Compute the determinant.
!
  if ( job/10 /= 0 ) then

     det(1) = 1.0_rkx
     det(2) = 0.0_rkx
     s = 10.0_rkx
     ii = 0

     do i = 1, n

        ii = ii + i

        det(1) = ap(ii)**2 * det(1)

        if ( det(1) == 0.0_rkx ) then
          exit
        end if

        do while ( det(1) < 1.0_rkx )
          det(1) = s * det(1)
          det(2) = det(2) - 1.0_rkx
        end do

        do while ( det(1) >= s )
          det(1) = det(1) / s
          det(2) = det(2) + 1.0_rkx
        end do

    end do

  end if
!
!  Compute inverse(R).
!
  if ( mod ( job, 10 ) /= 0 ) then

     kk = 0

     do k = 1, n

        k1 = kk + 1
        kk = kk + k
        ap(kk) = 1.0_rkx / ap(kk)
        t = - ap(kk)
        call sscal ( k-1, t, ap(k1), 1 )
        j1 = kk + 1
        kj = kk + k

        do j = k+1, n
           t = ap(kj)
           ap(kj) = 0.0_rkx
           call saxpy ( k, t, ap(k1), 1, ap(j1), 1 )
           j1 = j1 + j
           kj = kj + j
        end do

     end do
!
!  Form inverse(R) * transpose(inverse(R)).
!
     jj = 0

     do j = 1, n

        j1 = jj + 1
        jj = jj + j
        k1 = 1
        kj = j1

        do k = 1, j-1
           t = ap(kj)
           call saxpy ( k, t, ap(j1), 1, ap(k1), 1 )
           k1 = k1 + k
           kj = kj + 1
        end do

        t = ap(jj)
        call sscal ( j, t, ap(j1), 1 )

    end do

  end if

  return
end subroutine sppdi

subroutine sppfa ( ap, n, info )
!
!*******************************************************************************
!
!! SPPFA factors a real(rkx) symmetric positive definite matrix in packed form.
!
!
!  Discussion:
!
!    SPPFA is usually called by SPPCO, but it can be called
!    directly with a saving in time if RCOND is not needed.
!
!  Packed storage:
!
!    The following program segment will pack the upper
!    triangle of a symmetric matrix.
!
!      k = 0
!      do j = 1, n
!        do i = 1, j
!          k = k + 1
!          ap(k) = a(i,j)
!        end do
!      end do
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) AP(N*(N+1)/2).  On input, the packed form of a
!    symmetric matrix A.  The columns of the upper triangle are stored
!    sequentially in a one-dimensional array.  On output, an upper
!    triangular matrix R, stored in packed form, so that A = R'*R.
!
!    Input, integer N, the order of the matrix.
!
!    Output, integer INFO, error flag.
!    0, for normal return.
!    K, if the leading minor of order K is not positive definite.
!
  implicit none
!
  integer n
!
  real(rkx) ap((n*(n+1))/2)
  integer info
  integer j
  integer jj
  integer k
  integer kj
  integer kk
  real(rkx) s
  real(rkx) t
!
  info = 0
  jj = 0

  do j = 1, n

    s = 0.0_rkx
    kj = jj
    kk = 0

    do k = 1, j-1

      kj = kj + 1
      t = ap(kj) - sdot ( k-1, ap(kk+1), 1, ap(jj+1), 1 )
      kk = kk + k
      t = t / ap(kk)
      ap(kj) = t
      s = s + t * t

    end do

    jj = jj + j
    s = ap(jj) - s

    if ( s <= 0.0_rkx ) then
      info = j
      return
    end if

    ap(jj) = sqrt ( s )

  end do

  return
end subroutine sppfa

subroutine sppsl ( ap, n, b )
!
!*******************************************************************************
!
!! SPPSL solves a real(rkx) symmetric positive definite system factored by SPPCO or SPPFA.
!
!
!  Discussion:
!
!    To compute inverse(A) * C where C is a matrix with P columns
!
!      call sppco ( ap, n, rcond, z, info )
!
!      if ( rcond is too small .or. info /= 0 ) then
!        exit
!      end if
!
!      do j = 1, p
!        call sppsl ( ap, n, c(1,j) )
!      end do
!
!    A division by zero will occur if the input factor contains
!    a zero on the diagonal.  Technically this indicates
!    singularity but it is usually caused by improper subroutine
!    arguments.  It will not occur if the subroutines are called
!    correctly and INFO == 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) AP(N*(N+1)/2), the output from SPPCO or SPPFA.
!
!    Input, integer N, the order of the matrix.
!
!    Input/output, real(rkx) B(N).  On input, the right hand side.
!    On output, the solution.
!
  implicit none
!
  integer n
!
  real(rkx) ap((n*(n+1))/2)
  real(rkx) b(n)
  integer k
  integer kk
  real(rkx) t
!
  kk = 0

  do k = 1, n
    t = sdot ( k-1, ap(kk+1), 1, b(1), 1 )
    kk = kk + k
    b(k) = (b(k) - t)/ap(kk)
  end do

  do k = n, 1, -1
    b(k) = b(k) / ap(kk)
    kk = kk - k
    t = - b(k)
    call saxpy ( k-1, t, ap(kk+1), 1, b(1), 1 )
  end do

  return
end subroutine sppsl

subroutine sptsl ( n, d, e, b )
!
!*******************************************************************************
!
!! SPTSL solves a positive definite tridiagonal linear system.
!
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, integer N, the order of the matrix.
!
!    Input/output, real(rkx) D(N), on input the diagonal of the tridiagonal matrix.
!    On output, D is destroyed.
!
!    Input, real(rkx) E(N), the offdiagonal of the tridiagonal matrix in
!    entries E(1:N-1).
!
!    Input/output, real(rkx) B(N).  On input, the right hand side.
!    On output, the solution.
!
  implicit none
!
  integer n
!
  real(rkx) b(n)
  real(rkx) d(n)
  real(rkx) e(n)
  integer k
  integer kbm1
  integer ke
  integer kf
  integer kp1
  integer nm1d2
  real(rkx) t1
  real(rkx) t2
!
!  Check for 1 x 1 case.
!
  if ( n == 1 ) then
    b(1) = b(1) / d(1)
    return
  end if

  nm1d2 = ( n - 1 ) / 2

  if ( n > 2 ) then

    kbm1 = n - 1
!
!  Zero top half of subdiagonal and bottom half of superdiagonal.
!
    do k = 1, nm1d2
      t1 = e(k) / d(k)
      d(k+1) = d(k+1) - t1 * e(k)
      b(k+1) = b(k+1) - t1 * b(k)
      t2 = e(kbm1) / d(kbm1+1)
      d(kbm1) = d(kbm1) - t2 * e(kbm1)
      b(kbm1) = b(kbm1) - t2 * b(kbm1+1)
      kbm1 = kbm1 - 1
    end do

  end if

  kp1 = nm1d2 + 1
!
!  Clean up for possible 2 x 2 block at center.
!
  if ( mod ( n, 2 ) == 0 ) then
    t1 = e(kp1) / d(kp1)
    d(kp1+1) = d(kp1+1) - t1 * e(kp1)
    b(kp1+1) = b(kp1+1) - t1 * b(kp1)
    kp1 = kp1 + 1
  end if
!
!  Back solve starting at the center, going towards the top and bottom.
!
  b(kp1) = b(kp1)/d(kp1)

  if ( n > 2 ) then

    k = kp1 - 1
    ke = kp1 + nm1d2 - 1

    do kf = kp1, ke
      b(k) = ( b(k) - e(k) * b(k+1) ) / d(k)
      b(kf+1) = ( b(kf+1) - e(kf) * b(kf) ) / d(kf+1)
      k = k - 1
    end do

  end if

  if ( mod ( n, 2 ) == 0 ) then
    b(1) = ( b(1) - e(1) * b(2) ) / d(1)
  end if

  return
end subroutine sptsl

subroutine sqrdc ( a, lda, n, p, qraux, jpvt, work, job )
!
!*******************************************************************************
!
!! SQRDC computes the QR factorization of a real(rkx) rectangular matrix.
!
!
!  Discussion:
!
!    SQRDC uses Householder transformations.
!
!    Column pivoting based on the 2-norms of the reduced columns may be
!    performed at the user's option.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,P).  On input, the N by P matrix
!    whose decomposition is to be computed.  On output, A contains in
!    its upper triangle the upper triangular matrix R of the QR
!    factorization.  Below its diagonal A contains information from
!    which the orthogonal part of the decomposition can be recovered.
!    Note that if pivoting has been requested, the decomposition is not that
!    of the original matrix A but that of A with its columns permuted
!    as described by JPVT.
!
!    Input, integer LDA, the leading dimension of the array A.  LDA must
!    be at least N.
!
!    Input, integer N, the number of rows of the matrix A.
!
!    Input, integer P, the number of columns of the matrix A.
!
!    Output, real(rkx) QRAUX(P), contains further information required to recover
!    the orthogonal part of the decomposition.
!
!    Input/output, integer JPVT(P).  On input, JPVT contains integers that
!    control the selection of the pivot columns.  The K-th column A(*,K) of A
!    is placed in one of three classes according to the value of JPVT(K).
!      > 0, then A(K) is an initial column.
!      = 0, then A(K) is a free column.
!      < 0, then A(K) is a final column.
!    Before the decomposition is computed, initial columns are moved to
!    the beginning of the array A and final columns to the end.  Both
!    initial and final columns are frozen in place during the computation
!    and only free columns are moved.  At the K-th stage of the
!    reduction, if A(*,K) is occupied by a free column it is interchanged
!    with the free column of largest reduced norm.  JPVT is not referenced
!    if JOB == 0.  On output, JPVT(K) contains the index of the column of the
!    original matrix that has been interchanged into the K-th column, if
!    pivoting was requested.
!
!    Workspace, real(rkx) WORK(P).  WORK is not referenced if JOB == 0.
!
!    Input, integer JOB, initiates column pivoting.
!    0, no pivoting is done.
!    nonzero, pivoting is done.
!
  implicit none
!
  integer lda
  integer n
  integer p
!
  real(rkx) a(lda,p)
  integer jpvt(p)
  real(rkx) qraux(p)
  real(rkx) work(p)
  integer j
  integer job
  integer l
  integer lup
  integer maxj
  real(rkx) maxnrm
  real(rkx) nrmxl
  integer pl
  integer pu
  logical swapj
  real(rkx) t
  real(rkx) tt
!
  pl = 1
  pu = 0
!
!  If pivoting is requested, rearrange the columns.
!
  if ( job /= 0 ) then

    do j = 1, p

      swapj = jpvt(j) > 0

      if ( jpvt(j) < 0 ) then
        jpvt(j) = - j
      else
        jpvt(j) = j
      end if

      if ( swapj ) then

        if ( j /= pl ) then
          call sswap ( n, a(1,pl), 1, a(1,j), 1 )
        end if

        jpvt(j) = jpvt(pl)
        jpvt(pl) = j
        pl = pl + 1

      end if

    end do

    pu = p

    do j = p, 1, -1

      if ( jpvt(j) < 0 ) then

        jpvt(j) = - jpvt(j)

        if ( j /= pu ) then
          call sswap ( n, a(1,pu), 1, a(1,j), 1 )
          call i_swap ( jpvt(pu), jpvt(j) )
        end if

        pu = pu - 1

      end if

    end do

  end if
!
!  Compute the norms of the free columns.
!
  do j = pl, pu
    qraux(j) = snrm2 ( n, a(1,j), 1 )
  end do

  work(pl:pu) = qraux(pl:pu)
!
!  Perform the Householder reduction of A.
!
  lup = min ( n, p )

  do l = 1, lup
!
!  Bring the column of largest norm into the pivot position.
!
    if ( l >= pl .and. l < pu ) then

      maxnrm = 0.0_rkx
      maxj = l
      do j = l, pu
        if ( qraux(j) > maxnrm ) then
          maxnrm = qraux(j)
          maxj = j
        end if
      end do

      if ( maxj /= l ) then
        call sswap ( n, a(1,l), 1, a(1,maxj), 1 )
        qraux(maxj) = qraux(l)
        work(maxj) = work(l)
        call i_swap ( jpvt(maxj), jpvt(l) )
      end if

    end if
!
!  Compute the Householder transformation for column L.
!
    qraux(l) = 0.0_rkx

    if ( l /= n ) then

      nrmxl = snrm2 ( n-l+1, a(l,l), 1 )

      if ( nrmxl /= 0.0_rkx ) then

        if ( a(l,l) /= 0.0_rkx ) then
          nrmxl = sign ( nrmxl, a(l,l) )
        end if

        call sscal ( n-l+1, 1.0_rkx / nrmxl, a(l,l), 1 )
        a(l,l) = 1.0_rkx + a(l,l)
!
!  Apply the transformation to the remaining columns, updating the norms.
!
        do j = l + 1, p

          t = - sdot ( n-l+1, a(l,l), 1, a(l,j), 1 ) / a(l,l)
          call saxpy ( n-l+1, t, a(l,l), 1, a(l,j), 1 )

          if ( j >= pl .and. j <= pu ) then

            if ( qraux(j) /= 0.0_rkx ) then

              tt = 1.0_rkx - ( abs ( a(l,j) ) / qraux(j) )**2
              tt = max ( tt, 0.0_rkx )
              t = tt
              tt = 1.0_rkx + 0.05_rkx * tt * ( qraux(j) / work(j) )**2

              if ( tt /= 1.0_rkx ) then
                qraux(j) = qraux(j) * sqrt ( t )
              else
                qraux(j) = snrm2 ( n-l, a(l+1,j), 1 )
                work(j) = qraux(j)
              end if

            end if

          end if

        end do
!
!  Save the transformation.
!
        qraux(l) = a(l,l)
        a(l,l) = - nrmxl

      end if

    end if

  end do

  return
end subroutine sqrdc

subroutine sqrsl ( a, lda, n, k, qraux, y, qy, qty, b, rsd, ab, job, info )
!
!*******************************************************************************
!
!! SQRSL computes transformations, projections, and least squares solutions.
!
!
!  Discussion:
!
!    SQRSL requires the output of SQRDC.
!
!    For K <= min(N,P), let AK be the matrix
!
!      AK = ( A(JPVT(1)), A(JPVT(2)), ..., A(JPVT(K)) )
!
!    formed from columns JPVT(1), ..., JPVT(K) of the original
!    N by P matrix A that was input to SQRDC.  If no pivoting was
!    done, AK consists of the first K columns of A in their
!    original order.  SQRDC produces a factored orthogonal matrix Q
!    and an upper triangular matrix R such that
!
!      AK = Q * (R)
!               (0)
!
!    This information is contained in coded form in the arrays
!    A and QRAUX.
!
!    The parameters QY, QTY, B, RSD, and AB are not referenced
!    if their computation is not requested and in this case
!    can be replaced by dummy variables in the calling program.
!    To save storage, the user may in some cases use mod_the same
!    array for different parameters in the calling sequence.  A
!    frequently occuring example is when one wishes to compute
!    any of B, RSD, or AB and does not need Y or QTY.  In this
!    case one may identify Y, QTY, and one of B, RSD, or AB, while
!    providing separate arrays for anything else that is to be
!    computed.
!
!    Thus the calling sequence
!
!      call sqrsl ( a, lda, n, k, qraux, y, dum, y, b, y, dum, 110, info )
!
!    will result in the computation of B and RSD, with RSD
!    overwriting Y.  More generally, each item in the following
!    list contains groups of permissible identifications for
!    a single calling sequence.
!
!      1. (Y,QTY,B) (RSD) (AB) (QY)
!
!      2. (Y,QTY,RSD) (B) (AB) (QY)
!
!      3. (Y,QTY,AB) (B) (RSD) (QY)
!
!      4. (Y,QY) (QTY,B) (RSD) (AB)
!
!      5. (Y,QY) (QTY,RSD) (B) (AB)
!
!      6. (Y,QY) (QTY,AB) (B) (RSD)
!
!    In any group the value returned in the array allocated to
!    the group corresponds to the last member of the group.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) A(LDA,P), contains the output of SQRDC.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the number of rows of the matrix AK.  It must
!    have the same value as N in SQRDC.
!
!    Input, integer K, the number of columns of the matrix AK.  K
!    must not be greater than min(N,P), where P is the same as in the
!    calling sequence to SQRDC.
!
!    Input, real(rkx) QRAUX(P), the auxiliary output from SQRDC.
!
!    Input, real(rkx) Y(N), a vector to be manipulated by SQRSL.
!
!    Output, real(rkx) QY(N), contains Q * Y, if requested.
!
!    Output, real(rkx) QTY(N), contains Q' * Y, if requested.
!
!    Output, real(rkx) B(K), the solution of the least squares problem
!      minimize norm2 ( Y - AK * B),
!    if its computation has been requested.  Note that if pivoting was
!    requested in SQRDC, the J-th component of B will be associated with
!    column JPVT(J) of the original matrix A that was input into SQRDC.
!
!    Output, real(rkx) RSD(N), the least squares residual Y - AK * B,
!    if its computation has been requested.  RSD is also the orthogonal
!    projection of Y onto the orthogonal complement of the column space
!    of AK.
!
!    Output, real(rkx) AB(N), the least squares approximation Ak * B, if its
!    computation has been requested.  AB is also the orthogonal projection
!    of Y onto the column space of A.
!
!    Input, integer JOB, specifies what is to be computed.  JOB has
!    the decimal expansion ABCDE, with the following meaning:
!
!      if A /= 0, compute QY.
!      if B /= 0, compute QTY.
!      if C /= 0, compute QTY and B.
!      if D /= 0, compute QTY and RSD.
!      if E /= 0, compute QTY and AB.
!
!    Note that a request to compute B, RSD, or AB automatically triggers
!    the computation of QTY, for which an array must be provided in the
!    calling sequence.
!
!    Output, integer INFO, is zero unless the computation of B has
!    been requested and R is exactly singular.  In this case, INFO is the
!    index of the first zero diagonal element of R, and B is left unaltered.
!
  implicit none
!
  integer k
  integer lda
  integer n
!
  real(rkx) a(lda,*)
  real(rkx) ab(n)
  real(rkx) b(k)
  logical cab
  logical cb
  logical cqty
  logical cqy
  logical cr
  integer info
  integer j
  integer jj
  integer job
  integer ju
  real(rkx) qraux(*)
  real(rkx) qty(n)
  real(rkx) qy(n)
  real(rkx) rsd(n)
  real(rkx) t
  real(rkx) temp
  real(rkx) y(n)
!
!  set info flag.
!
  info = 0
!
!  Determine what is to be computed.
!
  cqy =        job / 10000         /= 0
  cqty = mod ( job,  10000 )       /= 0
  cb =   mod ( job,   1000 ) / 100 /= 0
  cr =   mod ( job,    100 ) /  10 /= 0
  cab =  mod ( job,     10 )       /= 0

  ju = min ( k, n-1 )
!
!  Special action when N = 1.
!
  if ( ju == 0 ) then

    if ( cqy ) then
      qy(1) = y(1)
    end if

    if ( cqty ) then
      qty(1) = y(1)
    end if

    if ( cab ) then
      ab(1) = y(1)
     end if

    if ( cb ) then

      if ( a(1,1) == 0.0_rkx ) then
        info = 1
      else
        b(1) = y(1) / a(1,1)
      end if

    end if

    if ( cr ) then
      rsd(1) = 0.0_rkx
    end if

    return

  end if
!
!  Set up to compute QY or QTY.
!
  if ( cqy ) then
    qy(1:n) = y(1:n)
  end if

  if ( cqty ) then
    qty(1:n) = y(1:n)
  end if
!
!  Compute QY.
!
  if ( cqy ) then

    do jj = 1, ju

      j = ju - jj + 1

      if ( qraux(j) /= 0.0_rkx ) then
        temp = a(j,j)
        a(j,j) = qraux(j)
        t = - sdot ( n-j+1, a(j,j), 1, qy(j), 1 ) / a(j,j)
        call saxpy ( n-j+1, t, a(j,j), 1, qy(j), 1 )
        a(j,j) = temp
      end if

    end do

  end if
!
!  Compute Q'*Y.
!
     if ( cqty ) then

        do j = 1, ju
           if ( qraux(j) /= 0.0_rkx ) then
              temp = a(j,j)
              a(j,j) = qraux(j)
              t = - sdot ( n-j+1, a(j,j), 1, qty(j), 1 ) / a(j,j)
              call saxpy ( n-j+1, t, a(j,j), 1, qty(j), 1 )
              a(j,j) = temp
           end if
        end do

     end if
!
!  Set up to compute B, RSD, or AB.
!
     if ( cb ) then
       b(1:k) = qty(1:k)
     end if

     if ( cab ) then
       ab(1:k) = qty(1:k)
     end if

     if ( cr .and. k < n ) then
       rsd(k+1:n) = qty(k+1:n)
     end if

     if ( cab .and. k+1 <= n ) then
       ab(k+1:n) = 0.0_rkx
     end if

     if ( cr ) then
       rsd(1:k) = 0.0_rkx
     end if
!
!  Compute B.
!
     if ( cb ) then

        do jj = 1, k

          j = k - jj + 1

          if ( a(j,j) == 0.0_rkx ) then
            info = j
            exit
          end if

          b(j) = b(j)/a(j,j)

          if ( j /= 1 ) then
            t = -b(j)
            call saxpy ( j-1, t, a(1,j), 1, b, 1 )
          end if

        end do

     end if

     if ( cr .or. cab ) then
!
!  Compute RSD or AB as required.
!
        do jj = 1, ju

           j = ju - jj + 1

           if ( qraux(j) /= 0.0_rkx ) then

              temp = a(j,j)
              a(j,j) = qraux(j)

              if ( cr ) then
                t = - sdot ( n-j+1, a(j,j), 1, rsd(j), 1 ) / a(j,j)
                call saxpy ( n-j+1, t, a(j,j), 1, rsd(j), 1 )
              end if

              if ( cab ) then
                t = - sdot ( n-j+1, a(j,j), 1, ab(j), 1 ) / a(j,j)
                call saxpy ( n-j+1, t, a(j,j), 1, ab(j), 1 )
              end if

              a(j,j) = temp

           end if

        end do

  end if

  return
end subroutine sqrsl

subroutine srot ( n, x, incx, y, incy, c, s )
!
!*******************************************************************************
!
!! SROT applies a plane rotation.
!
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input/output, real(rkx) X(*), one of the vectors to be rotated.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real(rkx) Y(*), one of the vectors to be rotated.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
!    Input, real(rkx) C, S, parameters (presumably the cosine and sine of
!    some angle) that define a plane rotation.
!
  implicit none
!
  real(rkx) c
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real(rkx) s
  real(rkx) stemp
  real(rkx) x(*)
  real(rkx) y(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    do i = 1, n
      stemp = c * x(i) + s * y(i)
      y(i) = c * y(i) - s * x(i)
      x(i) = stemp
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      stemp = c * x(ix) + s * y(iy)
      y(iy) = c * y(iy) - s * x(ix)
      x(ix) = stemp
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end subroutine srot

subroutine srotg ( sa, sb, c, s )
!
!*******************************************************************************
!
!! SROTG constructs a Givens plane rotation.
!
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input/output, real(rkx) SA, SB, ...
!
!    Output, real(rkx) C, S, ...
!
  implicit none
!
  real(rkx) c
  real(rkx) r
  real(rkx) roe
  real(rkx) s
  real(rkx) sa
  real(rkx) sb
  real(rkx) scale
  real(rkx) z
!
  if ( abs ( sa ) > abs ( sb ) ) then
    roe = sa
  else
    roe = sb
  end if

  scale = abs ( sa ) + abs ( sb )

  if ( scale == 0.0_rkx ) then
    c = 1.0_rkx
    s = 0.0_rkx
    r = 0.0_rkx
  else
    r = scale * sqrt ( ( sa / scale )**2 + ( sb / scale )**2 )
    r = sign ( 1.0_rkx, roe ) * r
    c = sa / r
    s = sb / r
  end if

  if ( abs ( c ) > 0.0_rkx .and. abs ( c ) <= s ) then
    z = 1.0_rkx / c
  else
    z = s
  end if

  sa = r
  sb = z

  return
end subroutine srotg

subroutine sscal ( n, sa, x, incx )
!
!*******************************************************************************
!
!! SSCAL scales a vector by a constant.
!
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real(rkx) SA, the multiplier.
!
!    Input/output, real(rkx) X(*), the vector to be scaled.
!
!    Input, integer INCX, the increment between successive entries of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer m
  integer n
  real(rkx) sa
  real(rkx) x(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 ) then

    m = mod ( n, 5 )

    x(1:m) = sa * x(1:m)

    do i = m+1, n, 5
      x(i)   = sa * x(i)
      x(i+1) = sa * x(i+1)
      x(i+2) = sa * x(i+2)
      x(i+3) = sa * x(i+3)
      x(i+4) = sa * x(i+4)
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    do i = 1, n
      x(ix) = sa * x(ix)
      ix = ix + incx
    end do

  end if

  return
end subroutine sscal

subroutine ssico ( a, lda, n, kpvt, rcond, z )
!
!*******************************************************************************
!
!! SSICO factors a real(rkx) symmetric matrix and estimates its condition.
!
!
!  Discussion:
!
!    If RCOND is not needed, SSIFA is slightly faster.
!
!    To solve A*X = B, follow SSICO by SSISL.
!
!    To compute inverse(A)*C, follow SSICO by SSISL.
!
!    To compute inverse(A), follow SSICO by SSIDI.
!
!    To compute determinant(A), follow SSICO by SSIDI.
!
!    To compute inertia(A), follow SSICO by SSIDI.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,N).  On input, the symmetric matrix to be
!    factored.  Only the diagonal and upper triangle are used.
!    On output, a block diagonal matrix and the multipliers which
!    were used to obtain it.  The factorization can be written A = U*D*U'
!    where U is a product of permutation and unit upper triangular
!    matrices, U' is the transpose of U, and D is block diagonal
!    with 1 by 1 and 2 by 2 blocks.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix.
!
!    Output, integer KPVT(N), pivot indices.
!
!    Output, real(rkx) RCOND, an estimate of the reciprocal condition of A.
!    For the system A*X = B, relative perturbations in A and B of size
!    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
!    If RCOND is so small that the logical expression
!      1.0_rkx + RCOND == 1.0_rkx
!    is true, then A may be singular to working precision.  In particular,
!    RCOND is zero if exact singularity is detected or the estimate underflows.
!
!    Output, real(rkx) Z(N), a work vector whose contents are usually unimportant.
!    If A is close to a singular matrix, then Z is an approximate null vector
!    in the sense that
!      norm(A*Z) = RCOND * norm(A) * norm(Z).
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  real(rkx) ak
  real(rkx) akm1
  real(rkx) anorm
  real(rkx) bk
  real(rkx) bkm1
  real(rkx) denom
  real(rkx) ek
  integer i
  integer info
  integer j
  integer k
  integer kp
  integer kps
  integer kpvt(n)
  integer ks
  real(rkx) rcond
  real(rkx) s
  real(rkx) ynorm
  real(rkx) z(n)
!
!  Find the norm of A, using only entries in the upper half of the matrix.
!
  do j = 1, n
    z(j) = sum ( abs ( a(1:j-1,j) ) )
    do i = 1, j-1
      z(i) = z(i) + abs ( a(i,j) )
    end do
  end do

  anorm = maxval ( z(1:n) )
!
!  Factor.
!
  call ssifa ( a, lda, n, kpvt, info )
!
!  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
!
!  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
!
!  The components of E are chosen to cause maximum local
!  growth in the elements of W where U*D*W = E.
!
!  The vectors are frequently rescaled to avoid overflow.
!
!  Solve U*D*W = E.
!
  ek = 1.0_rkx
  z(1:n) = 0.0_rkx
  k = n

  do while ( k /= 0 )

     if ( kpvt(k) < 0 ) then
       ks = 2
     else
       ks = 1
     end if

     kp = abs ( kpvt(k) )
     kps = k + 1 - ks

     if ( kp /= kps ) then
       call r_swap ( z(kps), z(kp) )
     end if

     if ( z(k) /= 0.0_rkx ) then
       ek = sign ( ek, z(k) )
     end if

     z(k) = z(k) + ek
     call saxpy ( k-ks, z(k), a(1,k), 1, z(1), 1 )

     if ( ks /= 1 ) then
       if ( z(k-1) /= 0.0_rkx ) then
         ek = sign ( ek, z(k-1) )
       end if
       z(k-1) = z(k-1) + ek
       call saxpy ( k-ks, z(k-1), a(1,k-1), 1, z(1), 1 )
     end if

     if ( ks /= 2 ) then

       if ( abs ( z(k) ) > abs ( a(k,k) ) ) then
         s = abs ( a(k,k) ) / abs ( z(k) )
         z(1:n) = s * z(1:n)
         ek = s * ek
       end if

       if ( a(k,k) /= 0.0_rkx ) then
         z(k) = z(k) / a(k,k)
       else
         z(k) = 1.0_rkx
       end if

     else

       ak = a(k,k) / a(k-1,k)
       akm1 = a(k-1,k-1) / a(k-1,k)
       bk = z(k) / a(k-1,k)
       bkm1 = z(k-1) / a(k-1,k)
       denom = ak * akm1 - 1.0_rkx
       z(k) = ( akm1 * bk - bkm1 ) / denom
       z(k-1) = ( ak * bkm1 - bk ) / denom

     end if

     k = k - ks

  end do

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )
!
!  Solve U'*Y = W.
!
  k = 1

  do while ( k <= n )

     ks = 1
     if ( kpvt(k) < 0 ) ks = 2

     if ( k /= 1 ) then

       z(k) = z(k) + sdot ( k-1, a(1,k), 1, z(1), 1 )

       if ( ks == 2 ) then
         z(k+1) = z(k+1) + sdot ( k-1, a(1,k+1), 1, z(1), 1 )
       end if

       kp = abs ( kpvt(k) )

       if ( kp /= k ) then
         call r_swap ( z(k), z(kp) )
       end if

     end if

     k = k + ks

  end do

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )

  ynorm = 1.0_rkx
!
!  Solve U*D*V = Y.
!
  k = n

  do while ( k /= 0 )

     if ( kpvt(k) < 0 ) then
       ks = 2
     else
       ks = 1
     end if

     if ( k /= ks ) then

       kp = abs ( kpvt(k) )
       kps = k + 1 - ks

       if ( kp /= kps ) then
         call r_swap ( z(kps), z(kp) )
       end if

       call saxpy ( k-ks, z(k), a(1,k), 1, z(1), 1 )

       if ( ks == 2 ) then
         call saxpy ( k-ks, z(k-1), a(1,k-1), 1, z(1), 1 )
       end if

     end if

     if ( ks /= 2 ) then

        if ( abs ( z(k) ) > abs ( a(k,k) ) ) then
          s = abs ( a(k,k) ) / abs ( z(k) )
          z(1:n) = s * z(1:n)
          ynorm = s * ynorm
        end if

        if ( a(k,k) /= 0.0_rkx ) then
          z(k) = z(k) / a(k,k)
        else
          z(k) = 1.0_rkx
        end if

     else

        ak = a(k,k) / a(k-1,k)
        akm1 = a(k-1,k-1) / a(k-1,k)
        bk = z(k) / a(k-1,k)
        bkm1 = z(k-1) / a(k-1,k)
        denom = ak * akm1 - 1.0_rkx
        z(k) = ( akm1 * bk - bkm1 ) / denom
        z(k-1) = ( ak * bkm1 - bk ) / denom

     end if

     k = k - ks

  end do

  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = s * ynorm
!
!  Solve U'*Z = V.
!
  k = 1

  do while ( k <= n )

     if ( kpvt(k) < 0 ) then
       ks = 2
     else
       ks = 1
     end if

     if ( k /= 1 ) then

        z(k) = z(k) + sdot ( k-1, a(1,k), 1, z(1), 1 )
        if ( ks == 2 ) then
          z(k+1) = z(k+1) + sdot ( k-1, a(1,k+1), 1, z(1), 1 )
        end if
        kp = abs ( kpvt(k) )

        if ( kp /= k ) then
          call r_swap ( z(k), z(kp) )
        end if

     end if

     k = k + ks

  end do
!
!  Make ZNORM = 1.0.
!
  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = s * ynorm

  if ( anorm /= 0.0_rkx ) then
    rcond = ynorm / anorm
  else
    rcond = 0.0_rkx
  end if

  return
end subroutine ssico

subroutine ssidi ( a, lda, n, kpvt, det, inert, work, job )
!
!*******************************************************************************
!
!! SSIDI computes the determinant, inertia and inverse of a real(rkx) symmetric matrix.
!
!
!  Discussion:
!
!    SSIDI uses the factors from SSIFA.
!
!    A division by zero may occur if the inverse is requested
!    and SSICO has set RCOND == 0.0_rkx or SSIFA has set INFO /= 0.
!
!    Variables not requested by JOB are not used.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,N).  On input, the output from SSIFA.
!    On output, the upper triangle of the inverse of the original matrix.
!    The strict lower triangle is never referenced.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer KPVT(N), the pivot vector from SSIFA.
!
!    Output, real(rkx) DET(2), the determinant of the original matrix.
!      determinant = DET(1) * 10.0**DET(2)
!    with 1.0_rkx <= abs ( DET(1) ) < 10.0_rkx or DET(1) = 0.0.
!
!    Output, integer INERT(3), the inertia of the original matrix.
!    INERT(1) = number of positive eigenvalues.
!    INERT(2) = number of negative eigenvalues.
!    INERT(3) = number of zero eigenvalues.
!
!    Workspace, real(rkx) WORK(N).
!
!    Input, integer JOB, specifies the tasks.
!    JOB has the decimal expansion ABC where
!    If C /= 0, the inverse is computed,
!    If B /= 0, the determinant is computed,
!    If A /= 0, the inertia is computed.
!    For example, JOB = 111 gives all three.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  real(rkx) ak
  real(rkx) akkp1
  real(rkx) akp1
  real(rkx) d
  real(rkx) det(2)
  logical dodet
  logical doert
  logical doinv
  integer inert(3)
  integer j
  integer jb
  integer job
  integer k
  integer kpvt(n)
  integer ks
  integer kstep
  real(rkx) t
  real(rkx), parameter :: ten = 10.0_rkx
  real(rkx) work(n)
!
  doinv = mod ( job,   10 )       /= 0
  dodet = mod ( job,  100 ) /  10 /= 0
  doert = mod ( job, 1000 ) / 100 /= 0

  if ( dodet .or. doert ) then

    if ( doert ) then
      inert(1) = 0
      inert(2) = 0
      inert(3) = 0
    end if

    if ( dodet ) then
      det(1) = 1.0_rkx
      det(2) = 0.0_rkx
    end if

    t = 0.0_rkx

    do k = 1, n

      d = a(k,k)
!
!  2 by 2 block.
!
!  use det (d  s)  =  (d/t * c - t) * t,  t = abs ( s )
!          (s  c)
!  to avoid underflow/overflow troubles.
!
!  Take two passes through scaling.  Use T for flag.
!
      if ( kpvt(k) <= 0 ) then

        if ( t == 0.0_rkx ) then
          t = abs ( a(k,k+1) )
          d = ( d / t ) * a(k+1,k+1) - t
        else
          d = t
          t = 0.0_rkx
        end if

      end if

      if ( doert ) then
        if ( d > 0.0_rkx ) then
          inert(1) = inert(1) + 1
        else if ( d < 0.0_rkx ) then
          inert(2) = inert(2) + 1
        else if ( d == 0.0_rkx ) then
          inert(3) = inert(3) + 1
        end if
      end if

      if ( dodet ) then

        det(1) = d * det(1)

        if ( det(1) /= 0.0_rkx ) then

          do while ( abs ( det(1) ) < 1.0_rkx )
            det(1) = ten * det(1)
            det(2) = det(2) - 1.0_rkx
          end do

          do while ( abs ( det(1) ) >= ten )
            det(1) = det(1) / ten
            det(2) = det(2) + 1.0_rkx
          end do

        end if

      end if

    end do

  end if
!
!  Compute inverse(A).
!
  if ( doinv ) then

     k = 1

     do while ( k <= n )

        if ( kpvt(k) >= 0 ) then
!
!  1 by 1.
!
           a(k,k) = 1.0_rkx / a(k,k)

           if ( k >= 2 ) then

              work(1:k-1) = a(1:k-1,k)

              do j = 1, k-1
                a(j,k) = sdot ( j, a(1,j), 1, work, 1 )
                call saxpy ( j-1, work(j), a(1,j), 1, a(1,k), 1 )
              end do

              a(k,k) = a(k,k) + sdot ( k-1, work, 1, a(1,k), 1 )

           end if

           kstep = 1
!
!  2 by 2.
!
        else

           t = abs ( a(k,k+1) )
           ak = a(k,k) / t
           akp1 = a(k+1,k+1) / t
           akkp1 = a(k,k+1) / t
           d = t * ( ak * akp1 - 1.0_rkx )
           a(k,k) = akp1 / d
           a(k+1,k+1) = ak / d
           a(k,k+1) = - akkp1 / d

           if ( k >= 2 ) then

              call scopy ( k-1, a(1,k+1), 1, work, 1 )

              do j = 1, k-1
                 a(j,k+1) = sdot ( j, a(1,j), 1, work, 1 )
                 call saxpy ( j-1, work(j), a(1,j), 1, a(1,k+1), 1 )
              end do

              a(k+1,k+1) = a(k+1,k+1) + sdot ( k-1, work, 1, a(1,k+1), 1 )
              a(k,k+1) = a(k,k+1) + sdot ( k-1, a(1,k), 1, a(1,k+1), 1 )
              call scopy ( k-1, a(1,k), 1, work, 1 )

              do j = 1, k-1
                 a(j,k) = sdot ( j, a(1,j), 1, work, 1 )
                 call saxpy ( j-1, work(j), a(1,j), 1, a(1,k), 1 )
              end do

              a(k,k) = a(k,k) + sdot ( k-1, work, 1, a(1,k), 1 )

           end if

           kstep = 2

        end if
!
!  Swap.
!
        ks = abs ( kpvt(k) )

        if ( ks /= k ) then

           call sswap ( ks, a(1,ks), 1, a(1,k), 1 )

           do jb = ks, k
             j = k + ks - jb
             call r_swap ( a(j,k), a(ks,j) )
           end do

           if ( kstep /= 1 ) then
             call r_swap ( a(ks,k+1), a(k,k+1) )
           end if

        end if

        k = k + kstep

    end do

  end if

  return
end subroutine ssidi

subroutine ssifa ( a, lda, n, kpvt, info )
!
!*******************************************************************************
!
!! SSIFA factors a real(rkx) symmetric matrix.
!
!
!  Discussion:
!
!     To solve A*X = B, follow SSIFA by SSISL.
!
!     To compute inverse(A)*C, follow SSIFA by SSISL.
!
!     To compute determinant(A), follow SSIFA by SSIDI.
!
!     To compute inertia(A), follow SSIFA by SSIDI.
!
!     To compute inverse(A), follow SSIFA by SSIDI.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) A(LDA,N).  On input, the symmetric matrix to be
!    factored.  Only the diagonal and upper triangle are used.
!    On output, a block diagonal matrix and the multipliers which
!    were used to obtain it.  The factorization can be written A = U*D*U'
!    where U is a product of permutation and unit upper triangular
!    matrices, U' is the transpose of U, and D is block diagonal
!    with 1 by 1 and 2 by 2 blocks.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix.
!
!    Output, integer KPVT(N), the pivot indices.
!
!    Output, integer INFO, error flag.
!    0, normal value.
!    K, if the K-th pivot block is singular.  This is not an error
!    condition for this subroutine, but it does indicate that SSISL
!    or SSIDI may divide by zero if called.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  real(rkx) absakk
  real(rkx) ak
  real(rkx) akm1
  real(rkx) alpha
  real(rkx) bk
  real(rkx) bkm1
  real(rkx) colmax
  real(rkx) denom
  integer imax
  integer imaxp1
  integer info
  integer j
  integer jj
  integer jmax
  integer k
  integer kpvt(n)
  integer kstep
  real(rkx) mulk
  real(rkx) mulkm1
  real(rkx) rowmax
  logical swap
  real(rkx) t
!
!  ALPHA is used in choosing pivot block size.
!
  alpha = ( 1.0_rkx + sqrt ( 17.0_rkx ) ) / 8.0_rkx

  info = 0
!
!  Main loop on K, which goes from N to 1.
!
  k = n

  do while ( k > 0 )

    if ( k == 1 ) then
      kpvt(1) = 1
      if ( a(1,1) == 0.0_rkx ) then
        info = 1
      end if
      return
    end if
!
!  This section of code determines the kind of
!  elimination to be performed.  When it is completed,
!  KSTEP will be set to the size of the pivot block, and
!  SWAP will be set to .true. if an interchange is required.
!
     absakk = abs ( a(k,k) )
!
!  Determine the largest off-diagonal element in column K.
!
     imax = isamax ( k-1, a(1,k), 1 )
     colmax = abs ( a(imax,k) )

     if ( absakk >= alpha * colmax ) then
       kstep = 1
       swap = .false.
     else
!
!  Determine the largest off-diagonal element in row IMAX.
!
        rowmax = 0.0_rkx
        imaxp1 = imax + 1

        do j = imaxp1, k
          rowmax = max ( rowmax, abs ( a(imax,j) ) )
        end do

        if ( imax /= 1 ) then
          jmax = isamax ( imax-1, a(1,imax), 1 )
          rowmax = max ( rowmax, abs ( a(jmax,imax) ) )
        end if

        if ( abs ( a(imax,imax) ) >= alpha * rowmax ) then
          kstep = 1
          swap = .true.
        else if ( absakk >= alpha * colmax * ( colmax / rowmax ) ) then
          kstep = 1
          swap = .false.
        else
          kstep = 2
          swap = ( imax /= k-1 )
        end if

     end if
!
!  Column K is zero.
!  Set INFO and iterate the loop.
!
     if ( max ( absakk, colmax ) == 0.0_rkx ) then

        kpvt(k) = k
        info = k
!
!  1 x 1 pivot block.
!
!  Perform an interchange.
!
     else if ( kstep /= 2 ) then

        if ( swap ) then

           call sswap ( imax, a(1,imax), 1, a(1,k), 1 )

           do jj = imax, k
             j = k + imax - jj
             call r_swap ( a(j,k), a(imax,j) )
           end do

        end if
!
!  Perform the elimination.
!
        do jj = 1, k-1
          j = k - jj
          mulk = - a(j,k) / a(k,k)
          t = mulk
          call saxpy ( j, t, a(1,k), 1, a(1,j), 1 )
          a(j,k) = mulk
        end do
!
!  Set the pivot array.
!
        if ( swap ) then
          kpvt(k) = imax
        else
          kpvt(k) = k
        end if
!
!  2 x 2 pivot block.
!
!  Perform an interchange.
!
     else

        if ( swap ) then

           call sswap ( imax, a(1,imax), 1, a(1,k-1), 1 )

           do jj = imax, k-1
             j = k-1 + imax - jj
             call r_swap ( a(j,k-1), a(imax,j) )
           end do

           t = a(k-1,k)
           a(k-1,k) = a(imax,k)
           a(imax,k) = t

        end if
!
!  Perform the elimination.
!
        if ( k-2 /= 0 ) then

           ak = a(k,k) / a(k-1,k)
           akm1 = a(k-1,k-1) / a(k-1,k)
           denom = 1.0_rkx - ak * akm1

           do jj = 1, k-2

              j = k-1 - jj
              bk = a(j,k) / a(k-1,k)
              bkm1 = a(j,k-1) / a(k-1,k)
              mulk = ( akm1 * bk - bkm1 ) / denom
              mulkm1 = ( ak * bkm1 - bk ) / denom
              t = mulk
              call saxpy ( j, t, a(1,k), 1, a(1,j), 1 )
              t = mulkm1
              call saxpy ( j, t, a(1,k-1), 1, a(1,j), 1 )
              a(j,k) = mulk
              a(j,k-1) = mulkm1

           end do

        end if
!
!  Set the pivot array.
!
        if ( swap ) then
          kpvt(k) = -imax
        else
          kpvt(k) = 1 - k
        end if

        kpvt(k-1) = kpvt(k)

     end if

     k = k - kstep

  end do

  return
end subroutine ssifa

subroutine ssisl ( a, lda, n, kpvt, b )
!
!*******************************************************************************
!
!! SSISL solves a real(rkx) symmetric system factored by SSIFA.
!
!
!  Discussion:
!
!    To compute inverse(A) * C where C is a matrix with P columns
!
!      call ssifa ( a, lda, n, kpvt, info )
!
!      if ( info == 0 ) then
!        do j = 1, p
!          call ssisl ( a, lda, n, kpvt, c(1,j) )
!        end do
!      end if
!
!    A division by zero may occur if the inverse is requested
!    and SSICO has set RCOND == 0.0_rkx or SSIFA has set INFO /= 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) A(LDA,N), the output from SSIFA.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer KPVT(N), the pivot vector from SSIFA.
!
!    Input/output, real(rkx) B(N).  On input, the right hand side.
!    On output, the solution.
!
  implicit none
!
  integer lda
  integer n
!
  real(rkx) a(lda,n)
  real(rkx) ak
  real(rkx) akm1
  real(rkx) b(n)
  real(rkx) bk
  real(rkx) bkm1
  real(rkx) denom
  integer k
  integer kp
  integer kpvt(n)
!
!  Loop backward applying the transformations and D inverse to B.
!
  k = n

  do while ( k > 0 )

     if ( kpvt(k) >= 0 ) then
!
!  1 x 1 pivot block.
!
        if ( k /= 1 ) then

           kp = kpvt(k)
!
!  Interchange.
!
           if ( kp /= k ) then
             call r_swap ( b(k), b(kp) )
           end if
!
!  Apply the transformation.
!
           call saxpy ( k-1, b(k), a(1,k), 1, b(1), 1 )

        end if
!
!  Apply D inverse.
!
        b(k) = b(k)/a(k,k)
        k = k - 1

     else
!
!  2 x 2 pivot block.
!
        if ( k /= 2 ) then

           kp = abs (kpvt(k))
!
!  Interchange.
!
           if ( kp /= k-1 ) then
              call r_swap ( b(k-1), b(kp) )
           end if
!
!  Apply the transformation.
!
           call saxpy ( k-2, b(k), a(1,k), 1, b(1), 1 )
           call saxpy ( k-2, b(k-1), a(1,k-1), 1, b(1), 1 )

        end if
!
!  Apply D inverse.
!
        ak = a(k,k) / a(k-1,k)
        akm1 = a(k-1,k-1) / a(k-1,k)
        bk = b(k) / a(k-1,k)
        bkm1 = b(k-1) / a(k-1,k)
        denom = ak * akm1 - 1.0_rkx
        b(k) = ( akm1 * bk - bkm1 ) / denom
        b(k-1) = ( ak * bkm1 - bk ) / denom
        k = k - 2

    end if

  end do
!
!  Loop forward applying the transformations.
!
  k = 1

  do while ( k <= n )

     if ( kpvt(k) >= 0 ) then
!
!  1 x 1 pivot block.
!
        if ( k /= 1 ) then
!
!  Apply the transformation.
!
           b(k) = b(k) + sdot ( k-1, a(1,k), 1, b(1), 1 )
           kp = kpvt(k)
!
!  Interchange.
!
           if ( kp /= k ) then
             call r_swap ( b(k), b(kp) )
           end if

        end if

        k = k + 1

     else
!
!  2 x 2 pivot block.
!
        if ( k /= 1 ) then
!
!  Apply the transformation.
!
           b(k) = b(k) + sdot ( k-1, a(1,k), 1, b(1), 1 )
           b(k+1) = b(k+1) + sdot ( k-1, a(1,k+1), 1, b(1), 1 )
           kp = abs (kpvt(k))
!
!  Interchange.
!
           if ( kp /= k ) then
             call r_swap ( b(k), b(kp) )
           end if

        end if

        k = k + 2

    end if

  end do

  return
end subroutine ssisl

subroutine sspco ( ap, n, kpvt, rcond, z )
!
!*******************************************************************************
!
!! SSPCO factors a real(rkx) symmetric matrix stored in packed form.
!
!
!  Discussion:
!
!    SSPCO uses elimination with symmetric pivoting and estimates
!    the condition of the matrix.
!
!    If RCOND is not needed, SSPFA is slightly faster.
!
!    To solve A*X = B, follow SSPCO by SSPSL.
!
!    To compute inverse(A)*C, follow SSPCO by SSPSL.
!
!    To compute inverse(A), follow SSPCO by SSPDI.
!
!    To compute determinant(A), follow SSPCO by SSPDI.
!
!    To compute inertia(A), follow SSPCO by SSPDI.
!
!  Packed storage:
!
!    The following program segment will pack the upper triangle of a
!    symmetric matrix.
!
!      k = 0
!      do j = 1, n
!        do i = 1, j
!          k = k + 1
!          ap(k) = a(i,j)
!        end do
!      end do
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) AP(N*(N+1)/2).  On input, the packed form of a
!    symmetric matrix A.  The columns of the upper triangle are stored
!    sequentially in a one-dimensional array.  On output, a block diagonal
!    matrix and the multipliers which were used to obtain it, stored in
!    packed form.  The factorization can be written A = U*D*U'
!    where U is a product of permutation and unit upper triangular
!    matrices, U' is the transpose of U, and D is block diagonal
!    with 1 by 1 and 2 by 2 blocks.
!
!    Input, integer N, the order of the matrix.
!
!    Output, integer KPVT(N), the pivot indices.
!
!    Output, real(rkx) RCOND, an estimate of the reciprocal condition of A.
!    For the system A*X = B, relative perturbations in A and B of size
!    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
!    If RCOND is so small that the logical expression
!      1.0_rkx + RCOND == 1.0_rkx
!    is true, then A may be singular to working precision.  In particular,
!    RCOND is zero if exact singularity is detected or the estimate underflows.
!
!    Output, real(rkx) Z(N) a work vector whose contents are usually unimportant.
!    If A is close to a singular matrix, then Z is an approximate null
!    vector in the sense that
!      norm(A*Z) = RCOND * norm(A) * norm(Z).
!
  implicit none
!
  integer n
!
  real(rkx) ak
  real(rkx) akm1
  real(rkx) anorm
  real(rkx) ap((n*(n+1))/2)
  real(rkx) bk
  real(rkx) bkm1
  real(rkx) denom
  real(rkx) ek
  integer i
  integer ij
  integer ik
  integer ikm1
  integer ikp1
  integer info
  integer j
  integer j1
  integer k
  integer kk
  integer km1k
  integer km1km1
  integer kp
  integer kps
  integer kpvt(n)
  integer ks
  real(rkx) rcond
  real(rkx) s
  real(rkx) t
  real(rkx) ynorm
  real(rkx) z(n)
!
!  Find norm of A using only upper half.
!
  j1 = 1
  do j = 1, n
    z(j) = sum ( abs ( ap(j1:j1+j-1) ) )
    ij = j1
    j1 = j1 + j
    do i = 1, j-1
      z(i) = z(i) + abs ( ap(ij) )
      ij = ij + 1
    end do
  end do

  anorm = maxval ( z(1:n) )
!
!  Factor.
!
  call sspfa ( ap, n, kpvt, info )
!
!  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
!
!  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
!
!  The components of E are chosen to cause maximum local
!  growth in the elements of W where U*D*W = E.
!
!  The vectors are frequently rescaled to avoid overflow.
!
!  Solve U*D*W = E.
!
  ek = 1.0_rkx
  z(1:n) = 0.0_rkx

  k = n
  ik = (n*(n - 1))/2

  do while ( k /= 0 )

     kk = ik + k
     ikm1 = ik - (k - 1)
     ks = 1
     if ( kpvt(k) < 0 ) ks = 2
     kp = abs (kpvt(k))
     kps = k + 1 - ks

     if ( kp /= kps ) then
        t = z(kps)
        z(kps) = z(kp)
        z(kp) = t
     end if

     if ( z(k) /= 0.0_rkx ) then
       ek = sign ( ek, z(k) )
     end if

     z(k) = z(k) + ek
     call saxpy ( k-ks, z(k), ap(ik+1), 1, z(1), 1 )

     if ( ks /= 1 ) then
        if ( z(k-1) /= 0.0_rkx ) then
          ek = sign ( ek, z(k-1) )
        end if
        z(k-1) = z(k-1) + ek
        call saxpy ( k-ks, z(k-1), ap(ikm1+1), 1, z(1), 1 )
     end if

     if ( ks /= 2 ) then

        if ( abs ( z(k) ) > abs ( ap(kk) ) ) then
           s = abs ( ap(kk) ) / abs ( z(k) )
           z(1:n) = s * z(1:n)
           ek = s * ek
        end if

        if ( ap(kk) /= 0.0_rkx ) then
          z(k) = z(k) / ap(kk)
        else
          z(k) = 1.0_rkx
        end if

     else

        km1k = ik + k - 1
        km1km1 = ikm1 + k - 1
        ak = ap(kk) / ap(km1k)
        akm1 = ap(km1km1) / ap(km1k)
        bk = z(k) / ap(km1k)
        bkm1 = z(k-1) / ap(km1k)
        denom = ak * akm1 - 1.0_rkx
        z(k) = ( akm1 * bk - bkm1 ) / denom
        z(k-1) = ( ak * bkm1 - bk ) / denom

     end if

     k = k - ks
     ik = ik - k
     if ( ks == 2 ) then
       ik = ik - (k + 1)
     end if

  end do

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )
!
!  Solve U'*Y = W.
!
  k = 1
  ik = 0

  do while ( k <= n )

     if ( kpvt(k) < 0 ) then
       ks = 2
     else
       ks = 1
     end if

     if ( k /= 1 ) then

        z(k) = z(k) + sdot ( k-1, ap(ik+1), 1, z(1), 1 )
        ikp1 = ik + k

        if ( ks == 2 ) then
          z(k+1) = z(k+1) + sdot ( k-1, ap(ikp1+1), 1, z(1), 1 )
        end if

        kp = abs (kpvt(k))

        if ( kp /= k ) then
           t = z(k)
           z(k) = z(kp)
           z(kp) = t
        end if

     end if

     ik = ik + k
     if ( ks == 2 ) ik = ik + (k + 1)
     k = k + ks

  end do

  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = 1.0_rkx
!
!  Solve U*D*V = Y.
!
  k = n

  ik = n * ( n - 1 ) / 2

  do while ( k > 0 )

     kk = ik + k
     ikm1 = ik - (k - 1)

     if ( kpvt(k) < 0 ) then
       ks = 2
     else
       ks = 1
     end if

     if ( k /= ks ) then

        kp = abs ( kpvt(k) )
        kps = k + 1 - ks

        if ( kp /= kps ) then
           t = z(kps)
           z(kps) = z(kp)
           z(kp) = t
        end if

        call saxpy ( k-ks, z(k), ap(ik+1), 1, z(1), 1 )

        if ( ks == 2 ) then
          call saxpy ( k-ks, z(k-1), ap(ikm1+1), 1, z(1), 1 )
        end if

     end if

     if ( ks /= 2 ) then

        if ( abs ( z(k) ) > abs ( ap(kk) ) ) then
           s = abs ( ap(kk) ) / abs ( z(k) )
           z(1:n) = s * z(1:n)
           ynorm = s * ynorm
        end if

        if ( ap(kk) /= 0.0_rkx ) then
          z(k) = z(k) / ap(kk)
        else
          z(k) = 1.0_rkx
        end if

     else

        km1k = ik + k - 1
        km1km1 = ikm1 + k - 1
        ak = ap(kk) / ap(km1k)
        akm1 = ap(km1km1) / ap(km1k)
        bk = z(k) / ap(km1k)
        bkm1 = z(k-1) / ap(km1k)
        denom = ak * akm1 - 1.0_rkx
        z(k) = ( akm1 * bk - bkm1 ) / denom
        z(k-1) = ( ak * bkm1 - bk ) / denom

     end if

     k = k - ks
     ik = ik - k
     if ( ks == 2 ) ik = ik - (k + 1)

  end do

  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = s * ynorm
!
!  Solve U'*Z = V.
!
  k = 1
  ik = 0

  do while ( k <= n )

     ks = 1
     if ( kpvt(k) < 0 ) ks = 2

     if ( k /= 1 ) then

        z(k) = z(k) + sdot ( k-1, ap(ik+1), 1, z(1), 1 )
        ikp1 = ik + k

        if ( ks == 2 ) then
          z(k+1) = z(k+1) + sdot ( k-1, ap(ikp1+1), 1, z(1), 1 )
        end if

        kp = abs ( kpvt(k) )

        if ( kp /= k ) then
           t = z(k)
           z(k) = z(kp)
           z(kp) = t
        end if

     end if

     ik = ik + k
     if ( ks == 2 ) ik = ik + (k + 1)
     k = k + ks

  end do
!
!  Make ZNORM = 1.0.
!
  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = s * ynorm

  if ( anorm /= 0.0_rkx ) then
    rcond = ynorm / anorm
  else
    rcond = 0.0_rkx
  end if

  return
end subroutine sspco

subroutine sspdi ( ap, n, kpvt, det, inert, work, job )
!
!*******************************************************************************
!
!! SSPDI computes the determinant, inertia and inverse of a real(rkx) symmetric matrix.
!
!
!  Discussion:
!
!    SSPDI uses the factors from SSPFA, where the matrix is stored in
!    packed form.
!
!    A division by zero will occur if the inverse is requested
!    and SSPCO has set RCOND == 0.0_rkx or SSPFA has set INFO /= 0.
!
!    Variables not requested by JOB are not used.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) AP(N*(N+1)/2).  On input, the output from SSPFA.
!    On output, the upper triangle of the inverse of the original matrix,
!    stored in packed form.  The columns of the upper triangle are stored
!    sequentially in a one-dimensional array.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer KPVT(N), the pivot vector from SSPFA.
!
!    Output, real(rkx) DET(2), the determinant of the original matrix.
!      determinant = DET(1) * 10.0**DET(2)
!    with 1.0_rkx <= abs ( DET(1) ) < 10.0_rkx or DET(1) = 0.0.
!
!    Output, integer INERT(3), the inertia of the original matrix.
!    INERT(1) = number of positive eigenvalues.
!    INERT(2) = number of negative eigenvalues.
!    INERT(3) = number of zero eigenvalues.
!
!    Workspace, real(rkx) WORK(N).
!
!    Input, integer JOB, has the decimal expansion ABC where:
!      if A /= 0, the inertia is computed,
!      if B /= 0, the determinant is computed,
!      if C /= 0, the inverse is computed.
!    For example, JOB = 111  gives all three.
!
  implicit none
!
  integer n
!
  real(rkx) ak
  real(rkx) akkp1
  real(rkx) akp1
  real(rkx) ap((n*(n+1))/2)
  real(rkx) d
  real(rkx) det(2)
  logical dodet
  logical doert
  logical doinv
  integer ij
  integer ik
  integer ikp1
  integer iks
  integer inert(3)
  integer j
  integer jb
  integer jk
  integer jkp1
  integer job
  integer k
  integer kk
  integer kkp1
  integer km1
  integer kpvt(n)
  integer ks
  integer ksj
  integer kskp1
  integer kstep
  real(rkx) t
  real(rkx), parameter :: ten = 10.0_rkx
  real(rkx) work(n)
!
  doinv = mod ( job,   10 )       /= 0
  dodet = mod ( job,  100 ) /  10 /= 0
  doert = mod ( job, 1000 ) / 100 /= 0

  if ( dodet .or. doert ) then

     if ( doert ) then
       inert(1) = 0
       inert(2) = 0
       inert(3) = 0
     end if

     if ( dodet ) then
       det(1) = 1.0_rkx
       det(2) = 0.0_rkx
     end if

     t = 0.0_rkx
     ik = 0

     do k = 1, n

        kk = ik + k
        d = ap(kk)
!
!  2 by 2 block
!  use det (d  s)  =  (d/t * c - t) * t,  t = abs ( s )
!          (s  c)
!  to avoid underflow/overflow troubles.
!
!  Take two passes through scaling.  Use T for flag.
!
        if ( kpvt(k) <= 0 ) then

           if ( t == 0.0_rkx ) then
              ikp1 = ik + k
              kkp1 = ikp1 + k
              t = abs ( ap(kkp1) )
              d = ( d / t ) * ap(kkp1+1) - t
           else
              d = t
              t = 0.0_rkx
           end if

        end if

        if ( doert ) then
           if ( d > 0.0_rkx ) then
             inert(1) = inert(1) + 1
           else if ( d < 0.0_rkx ) then
             inert(2) = inert(2) + 1
           else if ( d == 0.0_rkx ) then
             inert(3) = inert(3) + 1
           end if
        end if

        if ( dodet ) then

           det(1) = d * det(1)

           if ( det(1) /= 0.0_rkx ) then

             do while ( abs ( det(1) ) < 1.0_rkx )
               det(1) = ten * det(1)
               det(2) = det(2) - 1.0_rkx
             end do

             do while ( abs ( det(1) ) >= ten )
               det(1) = det(1) / ten
               det(2) = det(2) + 1.0_rkx
             end do

          end if

        end if

        ik = ik + k

    end do

  end if
!
!  Compute inverse(A).
!
  if ( doinv ) then

     k = 1
     ik = 0

     do while ( k <= n )

        km1 = k - 1
        kk = ik + k
        ikp1 = ik + k
        kkp1 = ikp1 + k

        if ( kpvt(k) >= 0) then
!
!  1 by 1.
!
           ap(kk) = 1.0_rkx / ap(kk)

           if ( k-1 >= 1 ) then

              call scopy ( k-1, ap(ik+1), 1, work, 1 )
              ij = 0

              do j = 1, k-1
                 jk = ik + j
                 ap(jk) = sdot ( j, ap(ij+1), 1, work, 1 )
                 call saxpy ( j-1, work(j), ap(ij+1), 1, ap(ik+1), 1 )
                 ij = ij + j
              end do

              ap(kk) = ap(kk) + sdot ( k-1, work, 1, ap(ik+1), 1 )

           end if

           kstep = 1

        else
!
!  2 by 2.
!
           t = abs ( ap(kkp1) )
           ak = ap(kk) / t
           akp1 = ap(kkp1+1) / t
           akkp1 = ap(kkp1) / t
           d = t * ( ak * akp1 - 1.0_rkx )
           ap(kk) = akp1 / d
           ap(kkp1+1) = ak / d
           ap(kkp1) = -akkp1 / d

           if ( km1 >= 1 ) then

              call scopy ( km1, ap(ikp1+1), 1, work, 1 )
              ij = 0

              do j = 1, km1
                 jkp1 = ikp1 + j
                 ap(jkp1) = sdot ( j, ap(ij+1), 1, work, 1 )
                 call saxpy ( j-1, work(j), ap(ij+1), 1, ap(ikp1+1), 1 )
                 ij = ij + j
              end do

              ap(kkp1+1) = ap(kkp1+1) + sdot ( km1, work, 1, ap(ikp1+1), 1 )
              ap(kkp1) = ap(kkp1) + sdot ( km1, ap(ik+1), 1, ap(ikp1+1), 1 )
              call scopy ( km1, ap(ik+1), 1, work, 1 )
              ij = 0

              do j = 1, km1
                 jk = ik + j
                 ap(jk) = sdot ( j, ap(ij+1), 1, work, 1 )
                 call saxpy ( j-1, work(j), ap(ij+1), 1, ap(ik+1), 1 )
                 ij = ij + j
              end do

              ap(kk) = ap(kk) + sdot ( km1, work, 1, ap(ik+1), 1 )

           end if

           kstep = 2

        end if
!
!  Swap.
!
        ks = abs ( kpvt(k) )

        if ( ks /= k ) then

           iks = ( ks * ( ks - 1 ) ) / 2
           call sswap ( ks, ap(iks+1), 1, ap(ik+1), 1 )
           ksj = ik + ks

           do jb = ks, k
              j = k + ks - jb
              jk = ik + j
              call r_swap ( ap(jk), ap(ksj) )
              ksj = ksj - (j - 1)
           end do

           if ( kstep /= 1 ) then
              kskp1 = ikp1 + ks
              call r_swap ( ap(kskp1), ap(kkp1) )
           end if

        end if

        ik = ik + k
        if (kstep == 2) ik = ik + k + 1
        k = k + kstep

     end do

  end if

  return
end subroutine sspdi

subroutine sspfa ( ap, n, kpvt, info )
!
!*******************************************************************************
!
!! SSPFA factors a real(rkx) symmetric matrix stored in packed form.
!
!
!  Discussion:
!
!    To solve A*X = B, follow SSPFA by SSPSL.
!
!    To compute inverse(A)*C, follow SSPFA by SSPSL.
!
!    To compute determinant(A), follow SSPFA by SSPDI.
!
!    To compute inertia(A), follow SSPFA by SSPDI.
!
!    To compute inverse(A), follow SSPFA by SSPDI.
!
!  Packed storage:
!
!    The following program segment will pack the upper triangle of a
!    symmetric matrix.
!
!      k = 0
!      do j = 1, n
!        do i = 1, j
!          k = k + 1
!          ap(k)  = a(i,j)
!        end do
!      end do
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) AP(N*(N+1)/2).  On input, the packed form of a symmetric
!    matrix A.  The columns of the upper triangle are stored sequentially
!    in a one-dimensional array.  On output, a block diagonal matrix and
!    the multipliers which were used to obtain it stored in packed form.
!    The factorization can be written A = U*D*U' where U is a product of
!    permutation and unit upper triangular matrices, U' is the transpose
!    of U, and D is block diagonal with 1 by 1 and 2 by 2 blocks.
!
!    Input, integer N, the order of the matrix.
!
!    Output, integer KPVT(N), the pivot indices.
!
!    Output, integer INFO, error flag.
!    0, normal value.
!    K, if the K-th pivot block is singular.  This is not an error
!    condition for this subroutine, but it does indicate that SSPSL or
!    SSPDI may divide by zero if called.
!
  implicit none
!
  integer n
!
  real(rkx) absakk
  real(rkx) ak
  real(rkx) akm1
  real(rkx) alpha
  real(rkx) ap((n*(n+1))/2)
  real(rkx) bk
  real(rkx) bkm1
  real(rkx) colmax
  real(rkx) denom
  integer ij
  integer ik
  integer ikm1
  integer im
  integer imax
  integer imaxp1
  integer imim
  integer imj
  integer imk
  integer info
  integer j
  integer jj
  integer jk
  integer jkm1
  integer jmax
  integer jmim
  integer k
  integer kk
  integer km1
  integer km1k
  integer km1km1
  integer kpvt(n)
  integer kstep
  real(rkx) mulk
  real(rkx) mulkm1
  real(rkx) rowmax
  logical swap
  real(rkx) t
!
!  ALPHA is used in choosing pivot block size.
!
  alpha = ( 1.0_rkx + sqrt ( 17.0_rkx ) ) / 8.0_rkx

  info = 0
!
!  Main loop on K, which goes from N to 1.
!
  k = n
  ik = (n*(n - 1))/2

10 continue
!
!  Leave the loop if K=0 or K=1.
!
  if ( k == 0 ) then
    return
  end if

  if ( k == 1 ) then
    kpvt(1) = 1
    if ( ap(1) == 0.0_rkx ) then
      info = 1
    end if
    return
  end if
!
!  This section of code determines the kind of elimination to be performed.
!  When it is completed, KSTEP will be set to the size of the pivot block,
!  and SWAP will be set to .true. if an interchange is required.
!
  km1 = k - 1
  kk = ik + k
  absakk = abs ( ap(kk) )
!
!  Determine the largest off-diagonal element in column K.
!
  imax = isamax ( k-1, ap(ik+1), 1 )
  imk = ik + imax
  colmax = abs ( ap(imk) )

  if ( absakk >= alpha * colmax ) then

    kstep = 1
    swap = .false.
!
!  Determine the largest off-diagonal element in row IMAX.
!
  else

    rowmax = 0.0_rkx
    imaxp1 = imax + 1
    im = imax * ( imax - 1 ) / 2
    imj = im + 2 * imax
    do j = imaxp1, k
      rowmax = max ( rowmax, abs ( ap(imj) ) )
      imj = imj + j
    end do

    if ( imax /= 1 ) then
      jmax = isamax ( imax-1, ap(im+1), 1 )
      jmim = jmax + im
      rowmax = max ( rowmax, abs ( ap(jmim) ) )
    end if

    imim = imax + im

    if ( abs ( ap(imim) ) >= alpha * rowmax ) then
      kstep = 1
      swap = .true.
    else if ( absakk >= alpha * colmax * ( colmax / rowmax ) ) then
      kstep = 1
      swap = .false.
    else
      kstep = 2
      swap = imax /= km1
    end if

  end if
!
!  Column K is zero.  Set INFO and iterate the loop.
!
     if ( max ( absakk, colmax ) == 0.0_rkx ) then
       kpvt(k) = k
       info = k
       go to 190
     end if

     if ( kstep /= 2 ) then
!
!  1 x 1 pivot block.
!
        if ( swap ) then
!
!  Perform an interchange.
!
           call sswap ( imax, ap(im+1), 1, ap(ik+1), 1 )
           imj = ik + imax

           do jj = imax, k
              j = k + imax - jj
              jk = ik + j
              t = ap(jk)
              ap(jk) = ap(imj)
              ap(imj) = t
              imj = imj - (j - 1)
           end do

        end if
!
!  Perform the elimination.
!
        ij = ik - (k - 1)

        do jj = 1, km1
           j = k - jj
           jk = ik + j
           mulk = -ap(jk)/ap(kk)
           t = mulk
           call saxpy ( j, t, ap(ik+1), 1, ap(ij+1), 1 )
           ap(jk) = mulk
           ij = ij - (j - 1)
        end do
!
!  Set the pivot array.
!
        if ( swap ) then
          kpvt(k) = imax
        else
          kpvt(k) = k
        end if

     else
!
!  2 x 2 pivot block.
!
        km1k = ik + k - 1
        ikm1 = ik - (k - 1)
!
!  Perform an interchange.
!
        if ( swap ) then

           call sswap ( imax, ap(im+1), 1, ap(ikm1+1), 1 )
           imj = ikm1 + imax

           do jj = imax, km1
              j = km1 + imax - jj
              jkm1 = ikm1 + j
              t = ap(jkm1)
              ap(jkm1) = ap(imj)
              ap(imj) = t
              imj = imj - ( j - 1 )
           end do

           t = ap(km1k)
           ap(km1k) = ap(imk)
           ap(imk) = t

        end if
!
!  Perform the elimination.
!
        if ( k-2 /= 0 ) then

           ak = ap(kk) / ap(km1k)
           km1km1 = ikm1 + k - 1
           akm1 = ap(km1km1) / ap(km1k)
           denom = 1.0_rkx - ak * akm1
           ij = ik - ( k - 1 ) - ( k - 2 )

           do jj = 1, k-2

              j = km1 - jj
              jk = ik + j
              bk = ap(jk) / ap(km1k)
              jkm1 = ikm1 + j
              bkm1 = ap(jkm1) / ap(km1k)
              mulk = ( akm1 * bk - bkm1 ) / denom
              mulkm1 = ( ak * bkm1 - bk ) / denom
              t = mulk
              call saxpy ( j, t, ap(ik+1), 1, ap(ij+1), 1 )
              t = mulkm1
              call saxpy ( j, t, ap(ikm1+1), 1, ap(ij+1), 1 )
              ap(jk) = mulk
              ap(jkm1) = mulkm1
              ij = ij - (j - 1)
          end do

        end if
!
!  Set the pivot array.
!
        kpvt(k) = 1 - k
        if (swap) kpvt(k) = -imax
        kpvt(k-1) = kpvt(k)

     end if

190  continue

     ik = ik - (k - 1)
     if ( kstep == 2 ) then
       ik = ik - (k - 2)
     end if

     k = k - kstep

  go to 10
end subroutine sspfa

subroutine sspsl ( ap, n, kpvt, b )
!
!*******************************************************************************
!
!! SSPSL solves the real(rkx) symmetric system factored by SSPFA.
!
!
!  Discussion:
!
!    To compute inverse(A) * C where C is a matrix with P columns:
!
!      call sspfa ( ap, n, kpvt, info )
!
!      if ( info /= 0 ) go to ...
!
!      do j = 1, p
!        call sspsl ( ap, n, kpvt, c(1,j) )
!      end do
!
!    A division by zero may occur if SSPCO has set RCOND == 0.0_rkx
!    or SSPFA has set INFO /= 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) AP(N*(N+1)/2), the output from SSPFA.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer KPVT(N), the pivot vector from SSPFA.
!
!    Input/output, real(rkx) B(N).  On input, the right hand side.
!    On output, the solution.
!
  implicit none
!
  integer n
!
  real(rkx) ak
  real(rkx) akm1
  real(rkx) ap((n*(n+1))/2)
  real(rkx) b(n)
  real(rkx) bk
  real(rkx) bkm1
  real(rkx) denom
  integer ik
  integer ikm1
  integer ikp1
  integer k
  integer kk
  integer km1k
  integer km1km1
  integer kp
  integer kpvt(n)
!
!  Loop backward applying the transformations and D inverse to B.
!
  k = n
  ik = (n*(n - 1))/2

  do while ( k > 0 )

     kk = ik + k

     if ( kpvt(k) >= 0 ) then
!
!  1 x 1 pivot block.
!
        if (k /= 1) then

           kp = kpvt(k)
!
!  Interchange.
!
           if ( kp /= k ) then
             call r_swap ( b(k), b(kp) )
           end if
!
!  Apply the transformation.
!
           call saxpy ( k-1, b(k), ap(ik+1), 1, b(1), 1 )

        end if
!
!  Apply D inverse.
!
        b(k) = b(k) / ap(kk)
        k = k - 1
        ik = ik - k

     else
!
!  2 x 2 pivot block.
!
        ikm1 = ik - (k - 1)

        if ( k /= 2 ) then

           kp = abs ( kpvt(k) )
!
!  Interchange.
!
           if ( kp /= k-1 ) then
              call r_swap ( b(k-1), b(kp) )
           end if
!
!  Apply the transformation.
!
           call saxpy ( k-2, b(k), ap(ik+1), 1, b(1), 1 )
           call saxpy ( k-2, b(k-1), ap(ikm1+1), 1, b(1), 1 )

        end if
!
!  Apply D inverse.
!
        km1k = ik + k - 1
        kk = ik + k
        ak = ap(kk) / ap(km1k)
        km1km1 = ikm1 + k - 1
        akm1 = ap(km1km1) / ap(km1k)
        bk = b(k) / ap(km1k)
        bkm1 = b(k-1) / ap(km1k)
        denom = ak * akm1 - 1.0_rkx
        b(k) = ( akm1 * bk - bkm1 ) / denom
        b(k-1) = ( ak * bkm1 - bk ) / denom
        k = k - 2
        ik = ik - (k + 1) - k

    end if

  end do
!
!  Loop forward applying the transformations.
!
  k = 1
  ik = 0

  do while ( k <= n )

     if ( kpvt(k) >= 0 ) then
!
!  1 x 1 pivot block.
!
        if ( k /= 1 ) then
!
!  Apply the transformation.
!
           b(k) = b(k) + sdot ( k-1, ap(ik+1), 1, b(1), 1 )
           kp = kpvt(k)
!
!  Interchange.
!
           if ( kp /= k ) then
             call r_swap ( b(k), b(kp) )
           end if

        end if

        ik = ik + k
        k = k + 1

     else
!
!  2 x 2 pivot block.
!
        if ( k /= 1 ) then
!
!  Apply the transformation.
!
           b(k) = b(k) + sdot ( k-1, ap(ik+1), 1, b(1), 1 )
           ikp1 = ik + k
           b(k+1) = b(k+1) + sdot ( k-1, ap(ikp1+1), 1, b(1), 1 )
           kp = abs ( kpvt(k) )
!
!  Interchange.
!
           if ( kp /= k ) then
             call r_swap ( b(k), b(kp) )
           end if

        end if

        ik = ik + k + k + 1
        k = k + 2

    end if

  end do

  return
end subroutine sspsl

subroutine ssvdc ( x, ldx, n, p, s, e, u, ldu, v, ldv, work, job, info )
!
!*******************************************************************************
!
!! SSVDC computes the singular value decomposition of a real(rkx) rectangular matrix.
!
!
!  Discussion:
!
!    SSVDC reduces a real(rkx) N by P matrix X to diagonal form by orthogonal
!    transformations U and V.  The diagonal elements S(I) are the singular
!    values of X.  The columns of U are the corresponding left singular
!    vectors, and the columns of V the right singular vectors.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) X(LDX,P).  On input, the matrix whose singular value
!    decomposition is to be computed.  On output, the matrix has been
!    destroyed.
!
!    Input, integer LDX, the leading dimension of the array X.
!    LDX must be at least N.
!
!    Input, integer N, the number of rows of the matrix.
!
!    Input, integer P, the number of columns of the matrix X.
!
!    Output, real(rkx) S(MM), where MM = min(N+1,P).  The first min(N,P) entries
!    of S contain the singular values of X arranged in descending
!    order of magnitude.
!
!    Output, real(rkx) E(P), ordinarily contains zeros.  However see the
!    discussion of INFO for exceptions.
!
!    Output, real(rkx) U(LDU,K).  If JOBA = 1 then K = N; if JOBA >= 2 then
!    K = min(N,P).  U contains the matrix of left singular vectors.
!    U is not referenced if JOBA = 0.  If N <= P or if JOBA = 2, then
!    U may be identified with X in the subroutine call.
!
!    Input, integer LDU, the leading dimension of the array U.
!    LDU must be at least N.
!
!    Output, real(rkx) V(LDV,P), the matrix of right singular vectors.
!    V is not referenced if JOB is 0.  If P <= N, then V may be identified
!    with X in the subroutine call.
!
!    Input, integer LDV, the leading dimension of the array V.
!    LDV must be at least P.
!
!    Workspace, real(rkx) WORK(N).
!
!    Input, integer JOB, controls the computation of the singular
!    vectors.  It has the decimal expansion AB with the following meaning:
!
!      A =  0, do not compute the left singular vectors.
!      A =  1, return the N left singular vectors in U.
!      A >= 2, return the first min(N,P) singular vectors in U.
!      B =  0, do not compute the right singular vectors.
!      B =  1, return the right singular vectors in V.
!
!    Output, integer INFO, status indicator.
!    The singular values (and their corresponding singular vectors)
!    S(INFO+1), S(INFO+2),...,S(M) are correct (here M=min(N,P)).
!    Thus if INFO is 0, all the singular values and their vectors are
!    correct.  In any event, the matrix B = U'*X*V is the bidiagonal matrix
!    with the elements of S on its diagonal and the elements of E on
!    its super-diagonal (U' is the transpose of U).  Thus the singular
!    values of X and B are the same.
!
  implicit none
!
  integer ldu
  integer ldv
  integer ldx
  integer n
  integer p
!
  real(rkx) b
  real(rkx) c
  real(rkx) cs
  real(rkx) e(p)
  real(rkx) el
  real(rkx) emm1
  real(rkx) f
  real(rkx) g
  integer info
  integer iter
  integer j
  integer job
  integer jobu
  integer k
  integer kase
  integer kk
  integer l
  integer ll
  integer lls
  integer ls
  integer lu
  integer m
  integer, parameter :: maxit = 30
  integer mm
  integer mm1
  integer nct
  integer nctp1
  integer ncu
  integer nrt
  integer nrtp1
  real(rkx) s(*)
  real(rkx) scale
  real(rkx) shift
  real(rkx) sl
  real(rkx) sm
  real(rkx) smm1
  real(rkx) sn
  real(rkx) t
  real(rkx) t1
  real(rkx) test
  real(rkx) u(ldu,*)
  real(rkx) v(ldv,p)
  logical wantu
  logical wantv
  real(rkx) work(n)
  real(rkx) x(ldx,p)
  real(rkx) ztest
!
!  Determine what is to be computed.
!
  wantu = .false.
  wantv = .false.
  jobu = mod ( job, 100 ) / 10

  if (jobu > 1) then
    ncu = min ( n, p )
  else
    ncu = n
  end if

  if ( jobu /= 0 ) then
    wantu = .true.
  end if

  if ( mod ( job, 10 ) /= 0 ) then
    wantv = .true.
  end if
!
!  Reduce X to bidiagonal form, storing the diagonal elements
!  in S and the super-diagonal elements in E.
!
  info = 0
  nct = min ( n-1, p )
  nrt = max ( 0, min ( p-2, n ) )
  lu = max ( nct, nrt )

  do l = 1, lu
!
!  Compute the transformation for the L-th column and
!  place the L-th diagonal in S(L).
!
     if ( l <= nct ) then

        s(l) = snrm2 ( n-l+1, x(l,l), 1 )

        if ( s(l) /= 0.0_rkx ) then
          if (x(l,l) /= 0.0_rkx ) then
            s(l) = sign ( s(l), x(l,l) )
          end if
          call sscal ( n-l+1, 1.0_rkx / s(l), x(l,l), 1 )
          x(l,l) = 1.0_rkx + x(l,l)
        end if

        s(l) = -s(l)

     end if

     do j = l+1, p
!
!  Apply the transformation.
!
        if ( l <= nct .and. s(l) /= 0.0_rkx ) then
           t = - sdot ( n-l+1, x(l,l), 1, x(l,j), 1 ) / x(l,l)
           call saxpy ( n-l+1, t, x(l,l), 1, x(l,j), 1 )
        end if
!
!  Place the L-th row of X into E for the
!  subsequent calculation of the row transformation.
!
        e(j) = x(l,j)
     end do
!
!  Place the transformation in U for subsequent back multiplication.
!
     if ( wantu .and. l <= nct ) then
        u(l:n,l) = x(l:n,l)
     end if

     if ( l <= nrt ) then
!
!  Compute the L-th row transformation and place the
!  L-th superdiagonal in E(L).
!
        e(l) = snrm2 ( p-l, e(l+1), 1 )

        if ( e(l) /= 0.0_rkx ) then
          if ( e(l+1) /= 0.0_rkx ) then
            e(l) = sign ( e(l), e(l+1) )
          end if
          call sscal ( p-l, 1.0_rkx / e(l), e(l+1), 1 )
          e(l+1) = 1.0_rkx + e(l+1)
        end if

        e(l) = -e(l)
!
!  Apply the transformation.
!
        if ( l+1 <= n .and. e(l) /= 0.0_rkx ) then

           work(l+1:n) = 0.0_rkx

           do j = l+1, p
              call saxpy ( n-l, e(j), x(l+1,j), 1, work(l+1), 1 )
           end do

           do j = l+1, p
              call saxpy ( n-l, -e(j)/e(l+1), work(l+1), 1, x(l+1,j), 1 )
           end do

        end if
!
!  Place the transformation in V for subsequent back multiplication.
!
        if ( wantv ) then
          v(l+1:p,l) = e(l+1:p)
        end if

     end if

  end do
!
!  Set up the final bidiagonal matrix of order M.
!
  m = min ( p, n+1 )
  nctp1 = nct + 1
  nrtp1 = nrt + 1

  if ( nct < p ) then
    s(nctp1) = x(nctp1,nctp1)
  end if

  if ( n < m ) then
    s(m) = 0.0_rkx
  end if

  if ( nrtp1 < m ) then
    e(nrtp1) = x(nrtp1,m)
  end if

  e(m) = 0.0_rkx
!
!  If required, generate U.
!
  if ( wantu ) then

    u(1:n,nctp1:ncu) = 0.0_rkx

    do j = nctp1, ncu
      u(j,j) = 1.0_rkx
    end do

    do ll = 1, nct

      l = nct - ll + 1

      if ( s(l) /= 0.0_rkx ) then

        do j = l+1, ncu
          t = - sdot ( n-l+1, u(l,l), 1, u(l,j), 1 ) / u(l,l)
          call saxpy ( n-l+1, t, u(l,l), 1, u(l,j), 1 )
        end do

        call sscal ( n-l+1, -1.0_rkx, u(l,l), 1 )
        u(l,l) = 1.0_rkx + u(l,l)
        u(1:l-1,l) = 0.0_rkx

      else

        u(1:n,l) = 0.0_rkx
        u(l,l) = 1.0_rkx

      end if

    end do

  end if
!
!  If it is required, generate V.
!
  if ( wantv ) then

     do ll = 1, p

        l = p - ll + 1

        if ( l <= nrt .and. e(l) /= 0.0_rkx ) then

           do j = l+1, p
              t = - sdot ( p-l, v(l+1,l), 1, v(l+1,j), 1 ) / v(l+1,l)
              call saxpy ( p-l, t, v(l+1,l), 1, v(l+1,j), 1 )
           end do

        end if

        v(1:p,l) = 0.0_rkx
        v(l,l) = 1.0_rkx

     end do

  end if
!
!  Main iteration loop for the singular values.
!
  mm = m
  iter = 0

  do while ( m > 0 )
!
!  If too many iterations have been performed, set flag and return.
!
     if ( iter >= maxit ) then
        info = m
        return
     end if
!
!  This section of the program inspects for
!  negligible elements in the S and E arrays.
!
!  On completion the variables KASE and L are set as follows:
!
!  KASE = 1     if S(M) and E(L-1) are negligible and L<M
!  KASE = 2     if S(L) is negligible and L<M
!  KASE = 3     if E(L-1) is negligible, L<M, and
!               S(L), ..., s(M) are not negligible (QR step).
!  KASE = 4     if E(M-1) is negligible (convergence).
!
     do ll = 1, m

        l = m - ll

        if ( l == 0 ) then
          exit
        end if

        test = abs ( s(l) ) + abs ( s(l+1) )
        ztest = test + abs ( e(l) )

        if ( ztest == test ) then
           e(l) = 0.0_rkx
           exit
        end if

     end do

     if ( l == m - 1 ) then

        kase = 4

     else

        do lls = l+1, m+1

           ls = m - lls + l+1

           if ( ls == l ) then
             exit
           end if

           test = 0.0_rkx
           if ( ls /= m ) then
             test = test + abs ( e(ls) )
           end if

           if ( ls /= l + 1 ) then
             test = test + abs ( e(ls-1) )
           end if

           ztest = test + abs ( s(ls) )

           if ( ztest == test ) then
              s(ls) = 0.0_rkx
              exit
           end if

        end do

        if ( ls == l ) then
          kase = 3
        else if ( ls == m ) then
          kase = 1
        else
          kase = 2
          l = ls
        end if

     end if

     l = l + 1
!
!  Deflate negligible S(M).
!
     if ( kase == 1 ) then

        mm1 = m - 1
        f = e(m-1)
        e(m-1) = 0.0_rkx

        do kk = l, mm1

           k = mm1 - kk + l
           t1 = s(k)
           call srotg ( t1, f, cs, sn )
           s(k) = t1

           if ( k /= l ) then
             f = - sn * e(k-1)
             e(k-1) = cs * e(k-1)
           end if

           if ( wantv ) then
             call srot ( p, v(1,k), 1, v(1,m), 1, cs, sn )
           end if

        end do
!
!  Split at negligible S(L).
!
     else if ( kase == 2 ) then

        f = e(l-1)
        e(l-1) = 0.0_rkx

        do k = l, m

           t1 = s(k)
           call srotg ( t1, f, cs, sn )
           s(k) = t1
           f = - sn * e(k)
           e(k) = cs * e(k)
           if ( wantu ) then
             call srot ( n, u(1,k), 1, u(1,l-1), 1, cs, sn )
           end if

        end do
!
!  Perform one QR step.
!
     else if ( kase == 3 ) then
!
!  Calculate the shift.
!
        scale = max ( abs ( s(m) ), abs ( s(m-1) ), abs ( e(m-1) ), &
                      abs ( s(l) ), abs ( e(l) ) )

        sm = s(m) / scale
        smm1 = s(m-1) / scale
        emm1 = e(m-1) / scale
        sl = s(l) / scale
        el = e(l) / scale
        b = ( ( smm1 + sm ) * ( smm1 - sm ) + emm1**2 ) / 2.0_rkx
        c = ( sm * emm1 )**2
        shift = 0.0_rkx

        if ( b /= 0.0_rkx .or. c /= 0.0_rkx ) then
           shift = sqrt ( b**2 + c )
           if ( b < 0.0_rkx ) then
             shift = - shift
           end if
           shift = c / ( b + shift )
        end if

        f = ( sl + sm ) * ( sl - sm ) - shift
        g = sl * el
!
!  Chase zeros.
!
        mm1 = m - 1

        do k = l, mm1

           call srotg ( f, g, cs, sn )

           if (k /= l) then
             e(k-1) = f
           end if

           f = cs * s(k) + sn * e(k)
           e(k) = cs * e(k) - sn * s(k)
           g = sn * s(k+1)
           s(k+1) = cs * s(k+1)

           if ( wantv ) then
             call srot ( p, v(1,k), 1, v(1,k+1), 1, cs, sn )
           end if

           call srotg ( f, g, cs, sn )
           s(k) = f
           f = cs * e(k) + sn * s(k+1)
           s(k+1) = -sn * e(k) + cs * s(k+1)
           g = sn * e(k+1)
           e(k+1) = cs * e(k+1)

           if ( wantu .and. k < n ) then
             call srot ( n, u(1,k), 1, u(1,k+1), 1, cs, sn )
           end if

        end do

        e(m-1) = f
        iter = iter + 1
!
!  Convergence.
!
     else if ( kase == 4 ) then
!
!  Make the singular value nonnegative.
!
        if ( s(l) < 0.0_rkx ) then
           s(l) = -s(l)
           if ( wantv ) then
             call sscal ( p, -1.0_rkx, v(1,l), 1 )
           end if
        end if
!
!  Order the singular value.
!
590     continue

        if ( l /= mm ) then

          if ( s(l) < s(l+1) ) then

            t = s(l)
            s(l) = s(l+1)
            s(l+1) = t

            if ( wantv .and. l < p ) then
              call sswap ( p, v(1,l), 1, v(1,l+1), 1 )
            end if

            if ( wantu .and. l < n ) then
              call sswap ( n, u(1,l), 1, u(1,l+1), 1 )
            end if

            l = l + 1

            go to 590

          end if

        end if

        iter = 0
        m = m - 1

    end if

  end do

  return
end subroutine ssvdc

subroutine sswap ( n, x, incx, y, incy )
!
!*******************************************************************************
!
!! SSWAP interchanges two vectors.
!
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input/output, real(rkx) X(*), one of the vectors to swap.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real(rkx) Y(*), one of the vectors to swap.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer m
  integer n
  real(rkx) stemp
  real(rkx) x(*)
  real(rkx) y(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    m = mod ( n, 3 )

    do i = 1, m
      stemp = x(i)
      x(i) = y(i)
      y(i) = stemp
    end do

    do i = m+1, n, 3

      stemp = x(i)
      x(i) = y(i)
      y(i) = stemp

      stemp = x(i + 1)
      x(i + 1) = y(i + 1)
      y(i + 1) = stemp

      stemp = x(i + 2)
      x(i + 2) = y(i + 2)
      y(i + 2) = stemp

    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      stemp = x(ix)
      x(ix) = y(iy)
      y(iy) = stemp
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end subroutine sswap

subroutine strco ( t, ldt, n, rcond, z, job )
!
!*******************************************************************************
!
!! STRCO estimates the condition of a real(rkx) triangular matrix.
!
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) T(LDT,N), the triangular matrix.  The zero elements of
!    the matrix are not referenced, and the corresponding elements of the
!    array can be used to store other information.
!
!    Input, integer LDT, the leading dimension of the array T.
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer JOB, indicates the shape of T:
!    0, T is lower triangular.
!    nonzero, T is upper triangular.
!
!    Output, real(rkx) RCOND, an estimate of the reciprocal condition of T.
!    For the system T*X = B, relative perturbations in T and B of size
!    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
!    If RCOND is so small that the logical expression
!      1.0_rkx + RCOND == 1.0_rkx
!    is true, then T may be singular to working precision.  In particular,
!    RCOND is zero if exact singularity is detected or the estimate underflows.
!
!    Output, real(rkx) Z(N) a work vector whose contents are usually unimportant.
!    If T is close to a singular matrix, then Z is an approximate null vector
!    in the sense that
!      norm(A*Z) = RCOND * norm(A) * norm(Z).
!
  implicit none
!
  integer ldt
  integer n
!
  real(rkx) ek
  integer i1
  integer j
  integer j1
  integer j2
  integer job
  integer k
  integer kk
  integer l
  logical lower
  real(rkx) rcond
  real(rkx) s
  real(rkx) sm
  real(rkx) t(ldt,n)
  real(rkx) tnorm
  real(rkx) w
  real(rkx) wk
  real(rkx) wkm
  real(rkx) ynorm
  real(rkx) z(n)
!
  lower = job == 0
!
!  Compute the 1-norm of T.
!
  tnorm = 0.0_rkx

  do j = 1, n

    if ( lower ) then
      l = n + 1 - j
      i1 = j
    else
      l = j
      i1 = 1
    end if

    tnorm = max ( tnorm, sum ( abs ( t(i1:i1+l-1,j) ) ) )

  end do
!
!  RCOND = 1/(norm(T)*(estimate of norm(inverse(T)))).
!
!  Estimate = norm(Z)/norm(Y) where T*Z = Y and T'*Y = E.
!
!  T' is the transpose of T.
!
!  The components of E are chosen to cause maximum local
!  growth in the elements of Y.
!
!  The vectors are frequently rescaled to avoid overflow.
!
!  Solve T'*Y = E.
!
  ek = 1.0_rkx
  z(1:n) = 0.0_rkx

  do kk = 1, n

     if ( lower ) then
       k = n + 1 - kk
     else
       k = kk
     end if

     if ( z(k) /= 0.0_rkx ) then
       ek = sign ( ek, -z(k) )
     end if

     if ( abs ( ek-z(k) ) > abs ( t(k,k) ) ) then
       s = abs ( t(k,k) ) / abs ( ek - z(k) )
       z(1:n) = s * z(1:n)
       ek = s * ek
     end if

     wk = ek - z(k)
     wkm = -ek - z(k)
     s = abs ( wk )
     sm = abs ( wkm )

     if ( t(k,k) /= 0.0_rkx ) then
       wk = wk / t(k,k)
       wkm = wkm / t(k,k)
     else
       wk = 1.0_rkx
       wkm = 1.0_rkx
     end if

     if ( kk /= n ) then

        if ( lower ) then
          j1 = 1
          j2 = k - 1
        else
          j1 = k + 1
          j2 = n
        end if

        do j = j1, j2
          sm = sm + abs ( z(j) + wkm * t(k,j) )
          z(j) = z(j) + wk * t(k,j)
          s = s + abs ( z(j) )
        end do

        if ( s < sm ) then
          w = wkm - wk
          wk = wkm
          z(j1:j2) = z(j1:j2) + w * t(k,j1:j2)
        end if

     end if

     z(k) = wk

  end do

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) )

  ynorm = 1.0_rkx
!
!  Solve T*Z = Y.
!
  do kk = 1, n

     if ( lower ) then
       k = kk
     else
       k = n + 1 - kk
     end if

     if ( abs ( z(k) ) > abs ( t(k,k) ) ) then
       s = abs ( t(k,k) ) / abs ( z(k) )
       z(1:n) = s * z(1:n)
       ynorm = s * ynorm
     end if

     if ( t(k,k) /= 0.0_rkx ) then
       z(k) = z(k) / t(k,k)
     else
       z(k) = 1.0_rkx
     end if

     if ( lower ) then
       i1 = k + 1
     else
       i1 = 1
     end if

     if ( kk < n ) then
       w = -z(k)
       call saxpy ( n-kk, w, t(i1,k), 1, z(i1), 1 )
     end if

  end do
!
!  Make ZNORM = 1.0.
!
  s = 1.0_rkx / sum ( abs ( z(1:n) ) )
  z(1:n) = s * z(1:n)
  ynorm = s * ynorm

  if ( tnorm /= 0.0_rkx ) then
    rcond = ynorm / tnorm
  else
    rcond = 0.0_rkx
  end if

  return
end subroutine strco

subroutine strdi ( t, ldt, n, det, job, info )
!
!***********************************************************************
!
!! STRDI computes the determinant and inverse of a real(rkx) triangular matrix.
!
!
!  Modified:
!
!    07 March 2001
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real(rkx) T(LDT,N).
!    On input, T contains the triangular matrix. The zero elements of the
!    matrix are not referenced, and the corresponding elements of the array
!    can be used to store other information.
!    On output, T contains the inverse matrix, if it was requested.
!
!    Input, integer LDT, the leading dimension of T.
!
!    Input, integer N, the order of the matrix.
!
!    Output, real(rkx) DET(2), the determinant of the matrix, if requested.
!    The determinant = DET(1) * 10.0**DET(2), with 1.0 <= abs ( DET(1) ) < 10.0,
!    or DET(1) == 0.
!
!    Input, integer JOB, specifies the shape of T, and the task.
!    010, inverse of lower triangular matrix.
!    011, inverse of upper triangular matrix.
!    100, determinant only.
!    110, determinant and inverse of lower triangular.
!    111, determinant and inverse of upper triangular.
!
!    Output, integer INFO.
!    If the inverse was requested, then
!    0, if the system was nonsingular;
!    nonzero, if the system was singular.
!
  implicit none
!
  integer ldt
  integer n
!
  real(rkx) det(2)
  integer i
  integer info
  integer j
  integer job
  integer k
  real(rkx) t(ldt,n)
  real(rkx) temp
  real(rkx), parameter :: ten = 10.0_rkx
!
!  Determinant.
!
  if ( job / 100 /= 0 ) then

    det(1) = 1.0_rkx
    det(2) = 0.0_rkx

    do i = 1, n

      det(1) = t(i,i) * det(1)

      if ( det(1) == 0.0_rkx ) then
        exit
      end if

      do while ( abs ( det(1) ) < 1.0_rkx )
        det(1) = ten * det(1)
        det(2) = det(2) - 1.0_rkx
      end do

      do while ( abs ( det(1) ) >= ten )
        det(1) = det(1) / ten
        det(2) = det(2) + 1.0_rkx
      end do

    end do

  end if

  if ( mod ( job / 10, 10 ) == 0 ) then
    return
  end if
!
!  Inverse of an upper triangular matrix.
!
  if ( mod ( job,10 ) /= 0 ) then

    info = 0

    do k = 1, n

      if ( t(k,k) == 0.0_rkx ) then
        info = k
        exit
      end if

      t(k,k) = 1.0_rkx / t(k,k)
      temp = -t(k,k)
      call sscal ( k-1, temp, t(1,k), 1 )

      do j = k + 1, n
        temp = t(k,j)
        t(k,j) = 0.0_rkx
        call saxpy ( k, temp, t(1,k), 1, t(1,j), 1 )
      end do

    end do
!
!  Inverse of a lower triangular matrix.
!
  else

    info = 0

    do k = n, 1, -1

      if ( t(k,k) == 0.0_rkx ) then
        info = k
        exit
      end if

      t(k,k) = 1.0_rkx / t(k,k)
      temp = -t(k,k)

      if ( k /= n ) then
        call sscal ( n-k, temp, t(k+1,k), 1 )
      end if

      do j = 1, k-1
        temp = t(k,j)
        t(k,j) = 0.0_rkx
        call saxpy ( n-k+1, temp, t(k,k), 1, t(k,j), 1 )
      end do

    end do

  end if

  return
end subroutine strdi

subroutine strsl ( t, ldt, n, b, job, info )
!
!*******************************************************************************
!
!! STRSL solves triangular linear systems.
!
!
!  Discussion:
!
!    STRSL can solve T * X = B or T' * X = B where T is a triangular
!    matrix of order N.
!
!    Here T' denotes the transpose of the matrix T.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, real(rkx) T(LDT,N), the matrix of the system.  The zero elements of
!    the matrix are not referenced, and the corresponding elements of the
!    array can be used to store other information.
!
!    Input, integer LDT, the leading dimension of the array T.
!
!    Input, integer N, the order of the matrix.
!
!    Input/output, real(rkx) B(N).  On input, the right hand side.
!    On output, the solution.
!
!    Input, integer JOB, specifies what kind of system is to be solved:
!
!    00, solve T * X = B, T lower triangular,
!    01, solve T * X = B, T upper triangular,
!    10, solve T'* X = B, T lower triangular,
!    11, solve T'* X = B, T upper triangular.
!
!    Output, integer INFO, singularity indicator.
!    0, the system is nonsingular.
!    nonzero, the index of the first zero diagonal element of T.
!
  implicit none
!
  integer ldt
  integer n
!
  real(rkx) b(n)
  integer icase
  integer info
  integer j
  integer jj
  integer job
  real(rkx) t(ldt,n)
  real(rkx) temp
!
!  Check for zero diagonal elements.
!
  do j = 1, n
    if ( t(j,j) == 0.0_rkx ) then
      info = j
      return
    end if
  end do

  info = 0
!
!  Determine the task and go to it.
!
     if ( mod ( job, 10 ) == 0 ) then
       icase = 1
     else
       icase = 2
     end if

     if ( mod ( job, 100 ) / 10 /= 0 ) then
       icase = icase + 2
     end if
!
!  Solve T*X=B for T lower triangular.
!
  if ( icase == 1 ) then

    b(1) = b(1) / t(1,1)

    do j = 2, n
      temp = -b(j-1)
      call saxpy ( n-j+1, temp, t(j,j-1), 1, b(j), 1 )
      b(j) = b(j) / t(j,j)
    end do
!
!  Solve T*X=B for T upper triangular.
!
  else if ( icase == 2 ) then

    b(n) = b(n) / t(n,n)

    do jj = 2, n
      j = n - jj + 1
      temp = -b(j+1)
      call saxpy ( j, temp, t(1,j+1), 1, b(1), 1 )
      b(j) = b(j) / t(j,j)
    end do
!
!  Solve T'*X=B for T lower triangular.
!
  else if ( icase == 3 ) then

    b(n) = b(n) / t(n,n)

    do jj = 2, n
      j = n - jj + 1
      b(j) = b(j) - sdot ( jj-1, t(j+1,j), 1, b(j+1), 1 )
      b(j) = b(j) / t(j,j)
    end do
!
!  Solve T'*X=B for T upper triangular.
!
  else if ( icase == 4 ) then

    b(1) = b(1) / t(1,1)

    do j = 2, n
      b(j) = b(j) - sdot ( j-1, t(1,j), 1, b(1), 1 )
      b(j) = b(j) / t(j,j)
    end do

  end if

  return
end subroutine strsl

end module linpack
! vim: tabstop=8 expandtab shiftwidth=2 softtabstop=2
